Resin monoterpene defenses decline within three widespread species of pine ( Pinus ) along a 1530‐m elevational gradient

Abstract

The elevational gradient in plant defense (EGPD) hypothesis posits that natural enemy pressures increase alongside temperature across elevational climatic gradients, thereby selecting for enhanced defenses at lower elevations while leaving plants less defended at higher elevations. Phylogenetically constrained tests of this hypothesis are uncommon, with tests focused on defenses of mature trees in natural settings being exceedingly rare. In the absence of this information, predicting the spatiotemporal dynamics of forest pests that preferentially attack mature trees is rendered more difficult. Tree properties such as age, growth rate, and size have all been correlated with levels of tree defenses against insect pests. Thus, we sought to disentangle the interacting influences of these properties from possible elevational climatic effects on monoterpene concentrations, diversity, and intraspecific dissimilarity of constitutive resin within three widespread pines (Pinus contorta, Pinus ponderosa, and Pinus flexilis) across a 1530‐m elevational transect in the Rocky Mountains, Colorado, USA. Overall, we found some support for the EGPD hypothesis within all three of the pine species studied. However, while elevation and tree properties were strong predictors of the variation in monoterpene concentration and moderate predictors of variation in diversity and dissimilarity in P. contorta, they were a relatively weak set of predictors in P. ponderosa and P. flexilis. Significant, (negative) elevational effects were found for total monoterpene concentration in P. contorta and monoterpene diversity in P. ponderosa, while elevation had a significant negative effect on monoterpene concentration and diversity, and a positive effect on dissimilarity in P. flexilis. Nevertheless, tree age, size, and/or growth measures had significant influences on one or more monoterpene measures in all three of the pine species indicating a need to parse effects of these factors within studies of tree defenses across environmental gradients. Our results offer some support for the EGPD hypothesis while suggesting the presence of a variable hierarchy of drivers of pine defense and resistance across congeners in shared environments—a situation that presents challenges for predicting spatiotemporal patterns in secondary chemistry across species from a common set of potential drivers.

Introduction

Trees possess defense mechanisms that enable some individuals to avoid or withstand pest attack (Trapp and Croteau , Franceschi et al. , Mumm and Hilker ). Tree defense phenotypes, while influenced by contemporary environmental conditions, are genetically controlled and have been shaped by abiotic and biotic selection pressures occurring at local and regional scales (Rosner and Hannrup , Holeski et al. , Sampedro et al. , Moreira et al. , Westbrook et al. , Zinkgraf et al. , Harper et al. ), as well as continental scales (Pearse and Hipp , Moreira et al. , Carrillo‐Gavilán et al. , Wang et al. ). As a result, tree defense phenotypes can vary substantially within and among populations (Moreira et al. , Ferrenberg et al. , Vogan and Schoettle , Zinkgraf et al. , Harper et al. ). This spatial variation in tree defenses, particularly as it relates to landscape features, has important implications for pest outbreaks (Holdenrieder et al. , Ferrenberg ).

Spatial variation in tree defense phenotypes can result from many dynamic, interacting mechanisms. At the scale of individual trees, defenses can increase or decrease across stages of ontogeny and phenology (Boege and Marquis , Barton and Koricheva ). For example, in lodgepole pines (Pinus contorta Dougl. ex Loud.), older stem material has been reported to be less chemically defended than younger material (Goodsman et al. ). In contrast, chemical defenses of white spruce (Picea glauca Moench) needles increase with age (Quiring ). With respect to both ontogeny and phenology, variations in tree size and growth rate have been linked to changes in defenses and risk of mortality (Ireland et al. , Bigler ). For example, in recent bark beetle outbreaks in western North America, larger‐diameter limber pines (Pinus flexilis James) suffered higher densities of bark beetle attacks (attacks/m² bark surface) than smaller trees (Ferrenberg and Mitton ) and larger‐diameter P. contorta were attacked at higher frequencies than smaller trees in some locations (Björklund and Lindgren ). At the same time, slower growth rates correlate with reduced production of resin ducts (specialized cells where resins are produced and stored) in the stems of numerous pine species, which increases susceptibility to bark beetle attack and mortality (Kane and Kolb , Gaylord et al. , Ferrenberg et al. , Hood et al. , Moreira et al. , Bentz et al. ). Importantly, Moreira et al. () found that resin defenses in the juvenile seedling stage of 17 pine species clustered into two divergent defensive syndromes with slow‐growing species (e.g., P. contorta and Pinus ponderosa Lawson) investing more heavily in constitutive defenses, while fast‐growing species have greater inducible defenses regardless of phylogenetic relatedness. Collectively, these findings suggest a set of interacting influences on pine defenses, whereby an individual tree's defenses are influenced by age, size, and growth rate within the context of species‐level adaptive strategies.

Tree defenses are also hypothesized to have evolved in response to variable biotic pressures across climatic gradients (Moles et al. ). Climatic suitability gradients for insects and pathogens are thought to influence the intensity of selection pressure placed upon plants, shaping defenses as a function of time and biogeography (Schemske et al. , Salazar and Marquis , Anstett et al. ). This idea is encapsulated in the longstanding latitudinal gradient in plant defenses (LGPD) hypothesis (Moles et al. ) and the more recent, yet conceptually related elevational gradient in plant defenses (EGPD) hypothesis (Mitton and Ferrenberg ). Both hypotheses posit that trees from warmer regions at lower latitudes and elevations experience more consistent pest pressures than trees from cooler regions at higher latitudes and elevations where pest ranges and active seasons are more frequently restricted by low temperatures (Salazar and Marquis ). Thus, in forest ecosystems, these selection gradients should result in greater tree defenses as a function of warmer temperatures (Moles et al. , Mitton and Ferrenberg , Pearse and Hipp , Moreira et al. , Abdala‐Roberts et al. , Wang et al. ).

In general, the LGPD and EGPD hypotheses have received mixed support across a range of plants. This may be due in part to experimental designs that fail to control for phylogeny, variation in defense traits across tree species, and variable feeding and attack strategies of different natural enemies (Anstett et al. ). Unfortunately, whether for conspecifics or congeners, comparisons of conifer defenses that examine a substantial elevational gradient are exceedingly rare (Hengxiao et al. , Smith , Abdala‐Roberts et al. ). In addition, to our knowledge, there are no previous studies of defenses across elevation that constrain phylogenetics and parse the influences of tree age, size, and growth rate on defense patterns despite their noted influences on tree defense and resistance levels.

The possible interacting influences of tree growth rate, size, age, and climatic gradients inspired us to examine conifer defensive traits along a 1532‐m elevational transect in the Rocky Mountains, Colorado, USA, using three widespread and ecologically important pine species of western North America: P. ponderosa var. scopulorum (Engelm), P. contorta var. latifolia (Engelm), and P. flexilis (James). Like all conifers, these species use resin as a primary defense against pests (Franceschi et al. , Mumm and Hilker ), investing significantly in constitutive resin defenses, particularly in stem material, as a resistance mechanism against specialized bark and wood‐boring insects that can cause substantial mortality during outbreaks and droughts (Ferrenberg et al. ). Pine resin consists of a suite of secondary metabolites, primarily diterpene acids dissolved in sesquiterpenes and monoterpenes. Variation in tree monoterpene phenotypes has been found to influence host selection and reproductive success of insect pests, as well as tree survival during pest outbreaks (Seybold et al. , Davis and Hofstetter , Raffa et al. , Taft et al. , b). Thus, in this study, we focused on the monoterpene components of resin collected from xylem and phloem of the tree stem given their importance in pine defense against various pests, including bark beetles. We hypothesized that monoterpene concentrations and diversity would be positively related to tree growth rates, and would decline with increasing tree age, size, and elevation. However, given potential tradeoffs among defenses and growth under optimal resource conditions, we considered an alternate hypothesis that higher growth rates would lead to decreased defenses as predicted by plant defense theory (Stamp ). We also hypothesized that intraspecific variation in monoterpene chemistry (i.e., level of dissimilarity among individuals within a species) would increase with elevation since less selection pressure by natural enemies in cooler regions should reduce directional and/or viability selection pressures against less‐effective defense phenotypes (Hoekstra et al. ). Our a priori prediction with regard to overall defense levels was that younger trees at lower elevations, with relatively fast growth rates, would be the most defended within our sampling region.

Materials and Methods

We tested the influence of tree age, size, and growth on tree defenses along climatic gradients by analyzing monoterpene chemistry of the resin of Pinus ponderosa, Pinus contorta, and Pinus flexilis across an elevational transect covering a total range of 1532 m from the lowest to highest sampling sites. The highest sampling site was located at 3330 m above sea level (asl; 40.0375, −105.5739; −0.3°C mean annual air temperature) and the lowest sampling site at 1795 m asl (40.1113, −105.3076; 9.6°C mean annual air temperature). Using the PRISM database (PRISM Climate Group ), we verified a strong correlation (R² = 0.99) between elevation and mean annual temperature (MAT) along this transect, with an average change of ±0.64°C for every 100 m of change in elevation (Appendix S1: Fig. S1). The highest, western‐most site (located within the Niwot Ridge LTER, Boulder County, Colorado, USA) and the lowest, eastern‐most site (located on private property near the city of Boulder, Colorado, USA) were separated by approximately 20 km of planar distance. The majority of sampling sites for intermediate elevations were located within a 10 km wide transect (north to south) lying between the high and low points of the east‐to‐west transect. Sampling sites were restricted to low and moderately sloping ridge‐tops and benches with a south–southeast aspect in locations where access was permitted by private landowners, the Niwot Ridge LTER, the University of Colorado's Mountain Research Station, the U.S. Forest Service, and Boulder County Parks and Open Space. Previous studies along or near this elevation gradient have revealed a decreasing stand disturbance frequency as elevation increases (Schuster et al. ) leading to a higher proportion of trees >100 yr for all three species in higher‐elevation stands. However, selective harvesting in the late 19th and early 20th century has led to a more homogenous age distribution in P. contorta compared to other conifers (Knowles and Grant ).

The species used here have large, often overlapping geographic ranges, yet have elevation ranges partitioned by climate. For instance, while all three species can be found growing together between 2579 and 2870 m asl, P. ponderosa dominates at lower elevations, while P. contorta and P. flexilis are more common at higher elevations. We attempted to sample each pine species across the majority of their elevational distributions, with P. ponderosa sampled from 1795 to 2854 m asl, P. contorta sampled from 2347 to 3330 m asl, and P. flexilis sampled from 2579 to 3330 m asl. Trees of each species were sampled between 2011 and 2014 from five or more locations along the elevational transect, with a diverse range of tree sizes included from along randomly oriented transects within each sampling site. The goal was to collect data on defenses from across a range of tree ages and growth rates while also effectively assessing elevational/climatic influences. In addition to characterizations of constitutive monoterpene chemistry, we measured tree diameter at breast height (dbh, 1.4 m above the ground surface) and collected an increment core from the sampled trees with a 12 mm diameter borer. However, of trees that supplied sufficient resin for chemical characterization, 12 individuals of P. contorta and 20 of P. ponderosa were lost after resin sampling but before physical measurements were taken due to wildfires or harvest. Monoterpene samples from trees without increment cores were retained only for analyses that did not include tree age or growth rates. Increment cores at breast height were used to estimate tree age using annual ring counts and annual radial growth (mm/yr) from the most recent five years (measured from 9600 dpi scans in Image J, Schneider et al. ). Radial growth rates were then used to calculate each tree's basal area increment (BAI) growth for the 5‐yr xylem ring interval. Basal area increment in this study was calculated as the percentage of a tree's total cross‐sectional area (mm²) comprising the most recent five years of xylem growth. The choice to use only a 5‐yr assessment of radial growth and BAI was made because resin chemistry—which can vary over time—was presumed to be more influenced by recent rather than historical tree growth and responses to the environment.

Resin collection and chemistry

Monoterpene types and concentrations were measured from resin collected in 3‐dram glass vials placed into a 6.4 mm diameter bore hole drilled at an upward 45° angle and approximately 25 mm deep into each tree's xylem. This approach was intended to capture constitutive resins but could allow for resin induction to occur, but importantly, we did not introduce pathogens or hormones to bore holes. Vials were tightly inserted to reduce contact between resin and the atmosphere, left in place for 48 h, and then capped and stored at −30°C until analysis. For laboratory analysis, resin samples were allowed to warm to room temperature before an aliquot of each sample was transferred to a 3‐dram glass vial which was promptly capped and weighed. Each sample was then dissolved in methyl tert‐butyl ether (C₅H₁₂O) containing either 5 mmol/L isobutylbenzene (C₁₀H₁₄) or 5 mmol/L cyclohexylbenzene (C₁₂H₁₆) as an internal standard. A Hewlett‐Packard HP6890 Series GC/MS (Hewlett‐Packard Company, Wilmington, Delaware, USA)equipped with a Restek Rtx‐5Sil MS column (Restek Corporation, Bellefonte, Pennsylvania, USA) was used to quantify monoterpene amounts. After split injection (2 μL, 34:1 ratio, 220°C), the oven temperature increased from an initial value of 40°C at 10°C/min until 200°C was reached (5‐min hold). During this time, helium flow rate was maintained at 1.6 mL/min. Identification of peaks was accomplished by comparing retention times to those of authentic standards. Calibration curves for quantification for the following monoterpenes were generated using authentic (standard) α‐pinene (Fluka Analytical, MilliporeSigma, St. Louis, Missouri, USA), β‐pinene (Fluka), α‐phellandrene (TCI America, Portland, Oregon, USA), β‐phellandrene (supplied by Dr. Ken Keefover‐Ring, University of Wisconsin), δ³‐carene (MilliporeSigma, St. Louis, Missouri, USA), limonene (Fluka), myrcene (Acros Organics, Thermo Fisher Scientific, Mullica Hill, New Jersey, USA), sabinene (Indofine Chemical Company, Inc., Hillsborough, New Jersey, USA), and terpinolene (TCI).

Statistical analyses

Relationships between tree age and size (dbh) were examined for each tree species with linear regression. Diameter at breast height for each species was log‐transformed to meet the assumptions of normality prior to regression analysis. Monoterpene concentration (mg monoterpenes/mg resin) and diversity (α‐diversity) calculated as the Shannon diversity index (H′) based on the composition and abundance of monoterpenes in each sample were compared among the three species using Welch's ANOVA for multigroup comparisons and Welch's t tests for pairwise comparisons to accommodate heteroscedasticity of group variances. Differences in the mean pairwise dissimilarity of monoterpene diversity (β‐diversity) calculated as the Bray–Curtis dissimilarity index of each species were compared to mean differences obtained by calculating all possible values of the test statistic under rearrangements of the labels on the observed data points, otherwise known as a permutation resample test. This approach determines whether observed differences in mean monoterpene dissimilarity among species differ from random. We used this approach to avoid violating assumptions of sample independence associated with other common tests. We used a stepwise logistic regression model to determine what proportion of resin samples could be correctly categorized to their species of derivation via monoterpene measures, with the best‐fit logistic regression model determined via the lowest value of the Bayesian information criterion.

We examined hypothesized relationships of total monoterpene concentration, diversity, dissimilarity, and individual monoterpene concentrations with elevation, tree age, size, growth rate (radial growth), and growth increment (BAI) via generalized linear models (GLMs). Generalized linear models were completed in JMP Pro 13.0.0 () using the adaptive lasso estimation method, an approach that computes parameter estimates by penalizing a weighted sum of the absolute values of the regression coefficients (Zou ). Best‐fit models were validated using the corrected Akaike information criterion. Data distribution types differed among GLMs. Concentration and diversity for P. contorta were modeled as a normal distribution using untransformed data; for all other GLMs dealing with total monoterpene concentration, diversity, and dissimilarity, the data were “normalized” via Johnson's Su transformations to approximate a normal curve and then modeled as normal distributions. For GLMs of individual monoterpene concentrations within each species, data were modeled as zero‐inflated gamma distributions (after verifying that exponential distributions matched the data) since not all trees sampled yielded each of the identified monoterpenes.

We used structural equation modeling (SEM; Grace ) to parameterize and evaluate the fit of hypothetical, causal models linking the monoterpene concentration and monoterpene diversity of constitutive resin to factors found to be significant in some GLMs (i.e., tree age, size, growth rate, and elevation). Similar to GLMs, monoterpene measures were only weakly explained at best by tree properties and elevation for P. ponderosa and P. flexilis resulting in poorly supported SEMs, leading us to discard models for these species. Structural equation models were performed in R version 3.1.2 (R Core Development Team 2014) with the lavaan package (Rosseel ). We used maximum‐likelihood parameter estimates and a mean‐adjusted Satorra–Bentler χ² (Satorra and Bentler ) to fit and test model fit, respectively. Because data used in the model were comprised of measures with various units, we calculated standardized coefficients as indicators of the magnitudes of the various paths in our SEM (Appendix S2: R code).

Results

We attempted to collect samples of constitutive resin from 106 individual trees of Pinus contorta, 118 individuals of Pinus ponderosa, and 106 individuals of Pinus flexilis across the elevational gradient. Of these trees, 63 P. contorta (59%), 108 P. ponderosa (92%), and 106 P. flexilis (100%) yielded a sufficient quantity of resin for further analysis. For trees yielding samples, the mean age of P. contorta was 104 yr with a range of 23–186, the mean age of P. ponderosa was 134 yr with a range of 21–424, and the mean age of P. flexilis was 117 with a range of 29–375. In P. contorta, mean tree size, measured as dbh, was 20.6 cm with a range of 11.3–36.1; in P. ponderosa, mean dbh was 20.1 with a range of 10.1–43.8; and in P. flexilis, mean dbh was 16.9 with a range of 3.7–39.2. For trees that supplied sufficient resin for analyses, dbh was positively related to tree age in P. ponderosa (R² = 0.23, P < 0.0001) and P. flexilis (R² = 0.24, P < 0.0001), but negatively related to tree age for P. contorta (R² = 0.12, P = 0.0103), indicating a likely bias of relatively small, older P. contorta trees supplying resin in our study, while relatively large, older trees did not since trees do not decrease in diameter as they age in the absence of physical damage or decay (Appendix S3: Fig. S1).

Total monoterpene diversity and concentration

Regardless of tree properties (age, growth rate, or size) or sample elevation, resin monoterpene concentration (mg monoterpene/mg resin) was significantly greater in P. contorta (0.21 ± 0.02) and P. ponderosa (0.20 ± 0.00) than in P. flexilis (0.14 ± 0.00; Welch's F = 42.5, df = 2, P < 0.0001; Fig. a). Relative abundance of monoterpenes differed among all three species (Appendix S4: Fig. S1). Best‐fit logistic regression retained a subset of six monoterpene concentrations (in order of their log worth for categorizing species: β‐phellandrene, limonene, sabinene, β‐pinene, δ³‐carene, and terpinolene) and correctly categorized 98.6% of resin samples by their tree species or origin (generalized R² = 0.99, Χ² = 572.1, P < 0.0001; Fig. ). All of the P. contorta samples were correctly classified, while two P. ponderosa were misclassified as being P. flexilis and two P. flexilis were misclassified as P. ponderosa, leading to an area under the receiver‐operating characteristic curve = 1.000 for P. contorta and = 0.999 for both P. ponderosa and P. flexilis (a value of 1.0 indicates a perfect grouping test, while a value of 0.5 indicates a test that performs no better than random assignment).

The diversity of resin monoterpenes (calculated as Shannon's H′, a measure of within‐sample diversity) significantly differed among all three species (Welch's F = 66.7, df = 2, P < 0.0001; Fig. b), with the highest diversity found in P. flexilis (1.57 ± 0.03), followed by P. ponderosa (1.31 ± 0.03) and P. contorta (1.09 ± 0.04). Mean pairwise dissimilarity in monoterpene diversity (a measure of average differences in monoterpene composition among samples) was also significantly different among the three species (Welch's F = 50.8, df = 2, P < 0.0001; Fig. c), with the highest level of monoterpene dissimilarity found in P. contorta (0.44 ± 0.01), followed by P. flexilis (0.41 ± 0.01) and P. ponderosa (0.33 ± 0.01).

Generalized regressions revealed that the combination of elevation and tree properties (age, size, and growth rates) were reasonably strong predictors of variation in monoterpene concentration (R² = 0.75) of P. contorta (Fig. ). Elevation was not retained in GLMs for diversity and dissimilarity of monoterpenes in P. contorta; however, tree age, size, and growth increment (BAI) were moderately good predictors of variation in diversity (R² = 0.31) and tree size and growth measures explained a smaller amount of variation in dissimilarity (R² = 0.23; Table ). Generalized linear models indicated that tree properties and elevation explained none of the variation in monoterpene concentration and dissimilarity in P. ponderosa (R² = 0.00 for both) but elevation, age, size, and BAI could explain a small amount of variation in monoterpene diversity (R² = 0.17; Table , Fig. a). For P. flexilis, GLMs also explained only a small portion of variation in concentration (R² = 0.17), diversity (R² = 0.16), and dissimilarity (R² = 0.19; Table ). However, elevation was retained as a significant factor in all three P. flexilis models and was the only factor retained in the GLM for monoterpene diversity (Fig. b–d).

Generalized regression models of monoterpene concentration (mg monoterpenes/mg resin), diversity (Shannon H′), and dissimilarity (Bray–Curtis dissimilarity) relationships to tree age, size (dbh), growth rate, and elevation in Pinus contorta, Pinus ponderosa, and Pinus flexilis

Measure	Factor	P. contorta	P. ponderosa	P. flexilis
χ 2	P‐value	Effect	R 2	χ 2	P‐value	Effect	R 2	χ 2	P‐value	Effect	R 2
Concentration	Age	30.6	<0.001	+	0.75	Removed		0.00	Removed		0.17
Size	20.9	<0.001	−	Removed		4.2	0.041	+
BAI	Removed		Removed		9.3	0.002	+
Radial growth	6.1	0.014	+	Removed		8.3	0.004	−
Elevation	28.7	<0.001	−	Removed		9.0	0.003	−
Diversity	Age	9.6	0.002	+	0.31	1.7	0.007	+	0.17	Removed		0.16
Size	5.2	0.023	−	3.3	0.070	+	Removed
BAI	2.0	0.161	+	0.1	0.743	−	Removed
Radial growth	Removed		Removed		Removed
Elevation	Removed		7.4	0.007	−	14.7	<0.001	−
Dissimilarity	Age	Removed		0.23	Removed		0.00	2.5	0.117	+	0.19
Size	1.3	0.250	+	Removed		Removed
BAI	2.2	0.139	+	Removed		0.6	0.450	−
Radial growth	1.2	0.281	−	Removed		Removed
Elevation	Removed		Removed		9.6	0.002	+

Notes

BAI, basal area increment; dbh, diameter at breast height. R² are for full generalized regression models; χ² = Wald chi‐square. Factors retained in each model were determined with adaptive lasso estimation, and best‐fit models were determined by corrected Akaike information criterion.

Individual monoterpene concentrations

Generalized linear models were also used to test the relationships of concentrations of individual monoterpenes (as opposed to total monoterpene concentration described above) with elevation and tree properties. Overall, GLMs suggested a significant influence of elevation on six of the 10 individual monoterpenes measured in P. contorta (P < 0.05; Appendix S5: Table S1), including three of the four most abundant (Appendix S6: Fig. S1). Importantly, tree age was the most influential factor in the variation in five of the 10 individual monoterpenes found in P. contorta (Appendix S5: Table S1). Generalized linear models for the four most abundant monoterpenes in P. contorta resulted in an R² = 0.70 for β‐phellandrene (most abundant), R² = 0.16 for δ3‐carene (second most abundant), R² = 0.00 for β‐pinene (third most abundant), and R² = 0.46 for α‐pinene (fourth most abundant). Generalized linear model fits ranged from R² = 0.23 to 0.07 for the remaining less abundant monoterpenes in P. contorta (Appendix S5: Table S1). Similarly, for P. ponderosa, elevation was retained as a significant factor in five of the 10 GLMs for individual monoterpenes. Generalized linear models for concentrations of the four most abundant individual monoterpenes (δ³‐carene, β‐pinene, myrcene, and α‐pinene, respectively) explained only a small portion of total variation with an R² = 0.13 for δ³‐carene (where elevation was retained as a significant factor; Appendix S6: Fig. S1) and R² = 0.00 for the other three most abundant monoterpenes (Appendix S5: Table S1). Generalized linear model fits ranged from R² = 0.64 to 0.00 for the remaining less abundant monoterpenes in P. ponderosa (Appendix S5: Table S1). Within P. flexilis, elevation was retained as a significant factor in just three of the 10 GLMs for individual monoterpenes and not retained for any of the four most abundant monoterpenes (δ³‐carene, limonene, α‐pinene, and sabinene, respectively). Overall, GLMs were poor fits for explaining variation in both abundant and more rare monoterpenes in P. flexilis with an R² = 0.04 for δ³‐carene, 0.00 for limonene, 0.06 for α‐pinene, and 0.02 sabinene and R² ranging from 0.15 to 0.00 for the remaining less abundant monoterpenes (Appendix S5: Table S1).

Structural equation models

The strong‐to‐moderate relationship of monoterpene concentration and diversity in P. contorta with tree properties and elevation in GLMs led us to propose and test hypothetical causal models (i.e., SEMs) of the interactive effects of tree age, size, growth rate, and elevation on P. contorta monoterpene concentration and diversity. Unlike GLMs where collinear predictors can lead to inflated model fit values, SEMs allow for predictor terms to interact via multiple specified paths. In our case, both SEMs were considered to be well fit to the observed data with the SEM for monoterpene concentration having a Satorra–Bentler χ² = 0.6, df = 1, and P = 0.43 (a small χ² and a large P‐value indicate an SEM that is well fit to observed data) and the SEM for monoterpene diversity having a Satorra–Bentler χ² = 0.4, df = 1, and P = 0.53. The SEMs for P. contorta explained 76% (R² = 0.76) of the variation in monoterpene concentration (Fig. a) and 33% (R² = 0.33) of the variation in monoterpene diversity (Fig. b). The SEM for monoterpene concentration in P. contorta resin revealed that (1) elevation has a negative direct effect on monoterpene concentration, (2) tree age has the strongest direct effect overall and is positively related to total monoterpene concentration, and has two paths of indirect influence via tree size and growth rate, and (3) increasing tree growth rates and sizes have conflicting direct effects on monoterpene concentrations, with trees that have grown faster (greater radial growth rate) over the most recent five years having greater concentrations, while larger‐diameter trees have lower concentrations (Fig. a). The SEM for monoterpene diversity also revealed that tree age has the strongest direct effect and increasing tree size has a negative direct effect, but revealed no significant links of monoterpene diversity to elevation (Fig. b).

Discussion

Pests rank among the leading causes of tree mortality and forest disturbance, with their impacts expected to increase due to global change pressures (Dietze and Matthes , Haynes et al. , Anderegg et al. , Kozlov et al. , Jules et al. ). Spatiotemporal patterns in tree defenses are an important feature of the landscape that affects forest pest dynamics (Ferrenberg ). Indeed, range expansion by insects has been linked to greater susceptibility of host plants within novel vs. historical ranges (Desurmont et al. ). Also, studies ranging in focus from herbs to conifers suggest that a diverse set of secondary chemicals is important for resisting natural enemies (Keefover‐Ring and Linhart , Iason et al. , Roach et al. , Richards et al. ). At the same time, other tree properties can influence risk of pest attack and mortality. In addition, resource and climatic gradients that can influence tree properties have also been implicated as controls on defenses in some species (O'Neill et al. , Moreira et al. , Ferrenberg ). A greater understanding of the drivers of spatiotemporal patterns in tree defenses is necessary to improve predictions of forest pest dynamics (Holdenrieder et al. , Ferrenberg ).

We characterized resin defenses in mature trees of three congeners occurring within shared environments across a 1530‐m elevational gradient. This design, while not an explicit test of population structure like a common garden, was necessary for supporting our goal of disentangling potentially interacting influences of age, growth, size, and elevation (a strong correlation of MATs in our study system; Appendix S1: Fig. S1) on defenses of mature trees of Pinus contorta, Pinus ponderosa, and Pinus flexilis. Space‐for‐time substitutions along natural climatic gradients to assess the influence of climate on plant ecology are now a common approach for the study of mature conifers given logistical constraints imposed by their long lifespans (De Frenne et al. ).

We found that the combined influences of tree properties and elevation explained a substantial amount of the variation in monoterpene concentration and a moderate amount of variation in monoterpene diversity and dissimilarity within P. contorta, but explained less variation in P. ponderosa and P. flexilis (Fig. , Table ). Nevertheless, we found support for our hypothesis that, by some measures, monoterpene defenses were negatively related to elevation in all three species (Table ). Elevation was the strongest direct negative influence on monoterpene concentration in P. contorta and P. flexilis (Table ). Similarly, elevation had the strongest direct negative effect on monoterpene diversity in P. ponderosa and P. flexilis, but had no influence on monoterpene concentration of P. ponderosa (Table ).

Tree growth rate, BAI, and size had important influences on some aspects of monoterpene defenses, but were not always retained as significant factors in models and had conflicting effects on defenses across species (Table ). Thus, generalizing the influence of growth and size on monoterpene defenses across these pine species is not immediately possible. Meanwhile, tree age was revealed as the strongest positive influence on monoterpene concentration and diversity in P. contorta and as a significant positive effect on diversity (but not concentration) in P. ponderosa, and dissimilarity (but not concentration or diversity) in P. flexilis. These findings oppose our hypothesis that defenses would decrease with age within all three of these congeners. Taken collectively with the elevation effects described above, these results suggest that defense chemistry is negatively related to elevational climatic gradients in all three of the pine species studied here, but that overall monoterpene concentrations and diversity are determined by a variable hierarchy of mechanisms among congeners. This result—while adding support to the EGPD hypothesis—suggests great challenges for predicting spatiotemporal patterns in tree secondary chemistry from the same set of factors across species.

In general, support for hypothesized climatic gradients in defenses, where plants are predicted to have greater defenses in warm regions, is mixed and a majority of studies refute their existence (Moles et al. , b). However, phylogenetically constrained experiments—those performed within the same species or genus—lend more support for the LGPD in both deciduous and coniferous trees (O'Neill et al. , Pearse and Hipp , Moreira et al. , Wang et al. , Zadworny et al. ). The EGPD hypothesis has also received support from both community‐level studies (Pellissier et al. , Rasmann et al. , Descombes et al. ) and phylogenetically constrained studies suggesting greater plant defenses at lower rather than higher elevation (Hengxiao et al. , Smith , O'Neill et al. , Anderson et al. , Dostálek et al. , Pellissier et al. ). However, phylogenetically constrained studies of tree defenses completed across substantial elevational gradients in field settings are exceedingly rare, possibly limited to our present study and those by Hengxiao et al. (), Smith (), and Abdala‐Roberts et al. (). These four studies diverge in their findings. Specially, while results from Hengxiao et al. () support the EGPD, Smith () found no change in monoterpene defenses of P. ponderosa (var. ponderosa Lawson) across a 1430‐m elevational transect of California's Sierra Nevada and Abdala‐Roberts et al. () reported an increase in leaf phenolics with elevation in English oak (Quercus robur) along an 830‐m elevation transect in northwestern Spain. This latter study also revealed strong influences of soil properties on tree defenses (Abdala‐Roberts et al. ). While phylogenetically constrained studies of insect damage on conifers are also rare, available evidence also suggests important influences of soil properties on pest impacts (Kosunen et al. ). Conflicting results among studies could be explained by the influence of abiotic drivers beyond temperature (e.g., soil properties) or by influences of biotic properties (e.g., tree age, growth rate, and size as demonstrated here) on tree defenses. Regardless, additional research is needed to better parse the role of climatic conditions from other potential drivers of plant defenses.

A strong influence of tree age on the number of pests attacking trees and in the amount of damage they cause has been reported in numerous studies (Boege et al. ). A meta‐analysis by Barton and Koricheva () found that mammalian herbivores consumed more tissue from older trees than younger, but found no clear age preference by phytophagous insects possibly due to variation in plant defense strategies and/or insect feeding types. As noted above, we found that tree age was a significant factor in some aspects of overall monoterpene defenses (i.e., concentration, diversity, and dissimilarity) for all three species (Table , Fig. ; Appendix S6: Fig. S1). Likewise, concentrations of seven individual monoterpenes within P. contorta, three within P. flexilis, and two within P. ponderosa were significantly influenced by tree age (Appendix S5: Table S1). Previous studies have also found that chemical defenses—particularly terpenoid compounds in needles—increased with age in various conifer species (Barton and Koricheva , Boege et al. ). However, Goodsman et al. () found that monoterpene concentrations were greater in younger xylem tissues than in older tissues of individual P. contorta trees. In addition to variation in tree ages and local environmental pressures across study sites, a plausible explanation for these seemingly contradictory results of our study with that of Goodsman et al. () possibly stems from their focus on phloem defenses vs. our focus on the combination of phloem and xylem defenses. Indeed, dramatic variation in levels and types of defenses across plant structures is a common feature of nearly all plant populations (Moore et al. , Hahn and Maron ). At the same time, we found a significant negative effect of P. contorta tree size on monoterpene concentration and diversity, which initially seems to be in conflict with an age‐related increase in defenses since older trees are generally expected to be larger than younger trees within a given population assuming similar growth conditions over time. Surprisingly, we found that the relationship between tree age and size was negative for P. contorta trees in our study, while the relationship between tree age and size was positive for P. ponderosa and P. flexilis. This negative relationship between tree age and size for P. contorta is clearly not mechanistic as trees do not decrease in diameter as they age in the absence of physical damage or decay. It instead appears to indicate that resin sampling was more successful from smaller, older P. contorta trees in this study system than from larger, older trees. Nevertheless, this negative relationship between age and size within our sample group, as well as the negative relationship between age and growth (Fig. ) in P. contorta that supplied resins in our study, might further explain the increase in monoterpene concentration and diversity with age, since slow‐growing conifers tend to produce higher concentrations of constitutive chemical defenses (Moreira et al. ). Also, evidence indicates that conifers exhibiting slower growth rates in early years of their lifespans tend ultimately to be smaller in size but live longer (Black et al. , Rötheli et al. , Bigler ) and suffer less pest damage later in life than trees with rapid early growth rates (Ruel and Whitham ).

An important consideration in the interpretation of our results is the finding that 43% of P. contorta, compared to 11% of P. ponderosa and 0% of P. flexilis, failed to produce a sufficient quantity of resin for monoterpene extraction upon wounding in our experiment. Smith () reported similarly low levels of resin yield in several populations of P. contorta, suggesting weak responses to mechanical wounding in the absence of pathogens. Whether the failure to pitch resin upon mechanical wounding indicates susceptibility to natural enemies is not entirely clear based on conflicting outcomes across tree and pest species. For example, while low resin flow could indicate poor production, trees can also over‐produce defensive resins and chemicals in response to relatively minor pest damage, thereby inducing their own mortality (Hulcr and Dunn ). Correcting yield bias across tree ages and sizes in P. contorta through alternate resin collection techniques that ensure yield from all targeted trees is an important next step in the study of conifer defenses as related to tree age, growth, and size across elevation.

Similar to results from previous comparisons of monoterpene defenses in these species (Smith , West et al. ), we found that overall, monoterpene communities of these pine species had little overlap (Fig. ; Appendix S4: Fig. S1), with interspecific variation far exceeding intraspecific variation regardless of elevational position or tree properties. The large differences among species were illustrated by the ability of logistic regression (based on monoterpene compositions and concentrations) to be >98% accurate for grouping resin samples within the proper species. Additionally, we found that monoterpene concentration was highest on average in P. contorta (Fig. ), while monoterpene diversity was greatest in P. flexilis (Fig. ). Increasing levels of both of these measures have been linked to greater resistance to natural enemies (Keefover‐Ring and Linhart , Davis and Hofstetter , Raffa ), making it difficult to assess which species is most or least defended in general. However, our finding that P. flexilis (subgenus Strobus) had lower monoterpene concentrations than congeners of the subgenus Pinus joins a result by Bentz et al. () who found P. flexilis to be generally “less defended” than co‐occurring Great Basin bristlecone pines (Pinus longaeva) and foxtail pines (Pinus balfouriana) despite these three species being members of the same subgenus. These consistent, interspecific variations in chemical defenses of pine species make it important to consider how deep evolutionary pressures may have shaped contemporary defense phenotypes and reinforce the value of experiments that consider how climate shapes defenses within species vs. approaches that compare defenses across species in the absence of knowledge regarding the pressures that gave rise to modern defense phenotypes (Ferrenberg et al. ).

Intraspecific variation in monoterpene chemistry was also substantial over the relatively small spatial scale of our study. Specifically, monoterpene dissimilarity—measured as the mean pairwise difference in monoterpene communities among each tree of a species—was roughly twofold greater in P. contorta than in P. ponderosa and P. flexilis (Fig. c). In other words, P. ponderosa and P. flexilis displayed less tree‐to‐tree variation in monoterpene defenses than did P. contorta. Our goal in measuring intraspecific monoterpene dissimilarity was to test the hypothesis that variation in defense phenotypes would increase with elevation as a reflection of less selection pressure by natural enemies. While levels of overall dissimilarity in P. contorta and P. ponderosa were not related to elevation, in P. flexilis dissimilarity in monoterpene profiles increased with elevation and tree age, such that older trees and those at higher elevation were less similar to neighboring conspecifics than younger and lower‐elevation trees. Previous work has unveiled similar levels of intraspecific variation in resin chemistry of the species studied here, but most typically over larger spatial scales than in our study (Forrest , Pollack and Dancik , Lindström et al. , Zavarin et al. , Smith , Pureswaran et al. , Clark et al. ). However, substantial variation in secondary resin chemistry among P. ponderosa was previously shown to occur near the lower end of the same elevational gradient used here (Latta et al. ) and over relatively small spatial scales in other locations of the U.S. Rocky Mountains (Davis and Hofstetter ).

Given that congeners occurring in the same environment must deal with similar abiotic and biotic pressures, large interspecific variation in resin chemistry alongside notable intraspecific variation suggests continued influences of deep evolutionary pressures alongside contemporary pressures on tree secondary chemistry. Elucidating the relationship between compositions and concentrations of constitutive and induced secondary chemicals is challenged by conflicting selection pressures at large and fine scales (Talluto and Benkman ). Local adaptations and low gene flow have been previously demonstrated for both P. ponderosa and P. flexilis populations along the same elevational transect used in our study (Schuster and Mitton ). Such local adaptation is likely driven by numerous factors. Evidence indicates that pressure by different natural enemies or abiotic factors can select for diverse monoterpene profiles (Iason et al. ) and lead to divergent plant defense strategies (Carmona and Fornoni , Moreira et al. ). Meta‐analytical results also indicate that secondary chemicals can have both detrimental and beneficial effects on insect pests of woody plants leading to conflicting adaptive significance over spatiotemporal scales (Smilanich et al. ). Additionally, profiles of monoterpenes and other chemicals in various coniferous tree species have been associated with a range of interactions with mammalian herbivores, insect pollinators, and resistance to fungal pathogens (Snyder , Langenheim , Latta and Mitton , Litvak and Monson , Witzell and Martín ). Similarly, a decrease in monoterpene concentrations (as shown in P. contorta and P. flexilis) or diversity (as shown in P. ponderosa and P. flexilis) at higher elevations, as well as the presence of greater dissimilarity in defense phenotypes with increasing elevation as found in P. flexilis, could be linked to physiological responses to lower atmospheric pressure, higher UV irradiance, and other factors unrelated to pest pressures (Bentz et al. ). Future studies focused on defenses beyond monoterpenes are also warranted since evidence from short‐lived herbaceous plants suggests the presence of variation in not only concentrations or density of defenses, but also in types of defenses expressed across climatic gradients (Dostálek et al. , Pellissier et al. ). In the case of P. flexilis populations used in this study, for example, trees can be defended against bark beetles via mechanical effects of bark texture even in the absence of resin flow (Ferrenberg and Mitton ). Finally, a focus on temperature in the absence of consideration for precipitation regime is potentially at the heart of incongruence results among studies testing the LGPD and EGPD hypotheses in coniferous trees. Work by O'Neill et al. () found that both constitutive and inducible defenses of spruce trees (Picea spp.) decreased with elevation and latitude in western Canada, but increased with aridity. A greater understanding of the drivers of interspecific and intraspecific variation in defensive monoterpene chemistry of the species studied here and others will be gained from additional experimental and observational efforts to disentangle the hierarchy of interacting influences on conifer defense expression.

Acknowledgments

We thank Boulder County Parks and Open Space, the Indian Peaks Wilderness Alliance, the John W. Marr Fund, and the University of Colorado's Department of Ecology and Evolutionary Biology for support. Publication of this article was funded by the University of Colorado Boulder Libraries Open Access Fund. Our manuscript was greatly improved by the comments of two anonymous reviewers. The authors declare no conflicts of interest.

Data Accessibility

Data from this study are available from figshare: https://doi.org/10.6084/m9.figshare.5386666.