The Legacy Effect of Mountain Pine Beetle Outbreaks on Chemical and Anatomical Defenses of Surviving Lodgepole Pine Trees

1. Introduction

Periodic outbreaks of bark beetles (Coleoptera: Curculionidae, Scolytinae) have caused extensive tree mortality, significantly altering forest structure and function [1,2,3] and shifting affected areas from carbon sinks to carbon sources [4,5,6]. For instance, recent outbreaks of the mountain pine beetle (MPB, Dendroctonus ponderosae Hopkins) have devastated millions of hectares of mature lodgepole pine (Pinus contorta) forests in western North America [1,2]. Despite this widespread mortality, certain individual trees?pine and non-pine species? have survived these outbreaks [7,8,9]. These surviving trees are crucial for maintaining ecological resilience within bark beetle-affected forests and may serve as a genetic reservoir of traits that confer resistance to future outbreaks, thereby contributing to the overall health and diversity of forest ecosystems [10,11,12,13,14]. Given their importance, numerous studies have reported on the impact of bark beetle outbreaks on these surviving trees, focusing on their post-outbreak defenses (chemical or anatomical) [7,8,9] and post-outbreak growing conditions [10,11,12,13,14]. While several studies have separately examined the impact of bark beetle outbreaks on anatomical and chemical defense mechanisms of the surviving trees [e.g. 7–9,15], few have explored the complex relationship between these defense modalities in the post-outbreak stands. This study aims to retrospectively examine the relationship between the constitutive terpenes in the current year’s phloem of surviving pine trees and the historical resin duct characteristics of the same trees before and after the MPB outbreak in lodgepole pine stands affected by the outbreak.

Interactions between conifers and bark beetles represent well-studied systems, integrating molecular, chemical, and ecological perspectives to elucidate the outcomes of these interactions [16,17,18]. Bark beetles feed, mate, and lay eggs within the subcortical tissues (phloem/cambial zone) of mature trees, making them prominent agents of mortality in coniferous forests [19]. When confronted with bark beetle attacks, conifer trees employ a multifaceted defense strategy, involving both defence compounds and anatomical structures [20,21,22,23,24,25,26]. Resin, a mixture of terpenes, serve as the primary defense against bark beetles for conifers. Monoterpenes and diterpenes, in particular, have received more attention due to their biological effects on bark beetles and their fungal associates [19,27,28,29,30]. Pine phloem, the feeding site for bark beetles, contains large quantities of monoterpenes and diterpenes, which are expressed constitutively and induced [19,23,31]. The production of these terpenes is supported in part by non-structural carbohydrates (NSC), such as soluble sugars and starch [32,33,34,35].

Anatomical defences of conifers involve primarily the synthesis and storage of resin in specialized structures known as resin cells and resin ducts [20,22]. These structures are interconnected axially and radially. Resin ducts are distributed throughout the primary and secondary xylem as well as the cortical tissues. Resin accumulated within these structures is transported to the point of bark beetle attacks. Consequently, the morphology of xylem resin ducts is often associated with the quantity of resin flow in response to wounds in the phloem [37,38,39,40,41]. In some cases, resin ducts are linked to tree survival against bark beetles, with larger duct size and higher duct numbers showing stronger relationships with survival [8,9,42,43,44].

Our earlier investigations demonstrated distinct anatomical defence traits in surviving trees following MPB outbreaks compared to those killed by MPBs [8,9]. Consequently, the surviving trees contained higher concentrations of monoterpenes [7]. However, the interaction between these chemical and anatomical defences remains unknown within the lodgepole pine system. Therefore, our primary objective is to retrospectively investigate whether the constitutive concentrations of post-outbreak monoterpenes and diterpenes in the phloem of surviving pine trees correlate with the historical resin duct characteristics, radial growth, and tree growth traits (age and size) before and after MPB outbreaks in post-outbreak lodgepole pine stands. Additionally, we explored potential relationships between defence traits (chemical or anatomical) and NSCs in the phloem and any trade-offs between monoterpene and diterpene defences following the MPB outbreak. This study aims to shed light on the complex relationship between various defence modalities and tree characteristics amidst the MPB outbreak, providing valuable insights into the legacy impact of bark beetles on host trees.

2. Materials and Methods

2.1. Site Selection

In July 2020, we identified five lodgepole pine-dominated forests in Jasper National Park, Alberta (Canada). Lodgepole pine constituted 45%-59% of the total live and dead tree density (height > 1.3 m) in sites 1-3 and 100% in sites 4 and 5 (Table S1). Douglas-fir (Pseudotsuga menziesii), white spruce (Picea glauca), and aspen (Populus tremuloides) were the other tree species in these sites. All sites had experienced a recent MPB outbreak, as documented in historical records provided by Parks Canada. The outbreaks occurred in sites 1-3 in 2015 and in sites 4 and 5 in 2014. The sites were located at least 500 meters apart. We gathered data on the on-site slope using a clinometer and aspect using a compass. Tree density at a site was calculated using an angle gauge with a basal area factor of 4.6 metric (20 standard). Tree species, health status (live or killed by MPB) and diameter at breast height (DBH, 1.3 m height) were measured for all trees within the plot (Table S1). Trees killed by MPB were identified through characteristic pitch tubes, brood emergence holes on the tree stems, blue stain and oviposition (maternal) galleries under the tree bark [7]. Conversely, the healthy lodgepole pines exhibited no visible symptoms of insect or pathogen attacks.

2.2. Sample Collection

We selected 12 pairs of beetle-killed and healthy lodgepole pine trees in each site, with DBH ranging from 18 to 22 cm (n=60 trees per category). We obtained two phloem samples using a 1.9 cm dia. leather hollow punch from each of the four cardinal directions at DBH on healthy trees (8 samples per tree). The punch was cleaned using 70% ethanol before each sampling. The samples were wrapped in labelled aluminum foil and immediately placed into liquid nitrogen. The samples were stored on dry ice during transportation from the field to the laboratory and put in a -40 °C freezer until chemical analysis.

Furthermore, we obtained one 5 mm and one 12 mm dia. increment cores at 1 m above ground from the north and south aspects of each health tree. The 5 mm dia. core was drilled to the pith of the tree and used to determine its age. The 12 mm dia. cores contained at least 20 years of annual rings dating back from the sampling year and were used to determine the anatomical defence and growth characteristics of trees in the last 15 years. The cores were stored in plastic drinking straws. To determine the year of the outbreak in each site, we took a wedge from the south-facing aspect of each beetle-killed tree at DBH. Wedges contained at least 20 years of growth, counting back from the year of death.

2.3. Chemical Analysis

After removing the outer bark from phloem samples, we ground the samples in liquid N with a mini bead beater at 2,200 rpm for three 20-s intervals. The eight samples from each tree were pooled together before analysis. The samples from each tree were divided into three groups for the following chemical analyses: monoterpenes (100 mg), diterpene (30 mg), and NSCs (30 mg). All these analyses followed the protocols reported by Cale et al. [45].

2.3.1. Monoterpene Analysis

At room temperature, 0.5 mL of hexane with pentadecane (0.004%) (internal standard) was mixed with 100 mg of sample. The mixture was vortexed for 30 s, sonicated for 10 min, and then centrifuged for 15 min at 18,213 g at 4 °C. This process was repeated twice for each sample. The extract was transferred to a 2 mL GC (Gas Chromatograph) vial and stored at -40 °C until analysis. Extracts were analyzed with a GC-MS (GC 7890A, Mass Spectrometry (MS) 5975C, Agilent Tech., Santa Clara, CA, USA) fitted with an Agilent HP-Chiral-20B column (30 m x 0.25 mm ID x 0.25 µm film). The carrier gas was helium at a flow rate of 1.0 mL min^-1. The temperature began at 50 °C and increased to 72 °C at a rate of 75 °C min^-1, where it was held for 0.5 min, then increased at 30 °C min^-1 to 90 °C and held for 2 min. It was then increased to 95 °C at 3 °C min^-1 and held for 1 min, then increased to 100 °C at 5 °C min^-1 and held for 0.5 min, then to 150 °C at 8 °C min^-1, and finally to 175 °C at 15 °C min^-1 and held for 1.5 min. The total run time for each sample was 23.6 min. We identified the hexane-extractable compounds, mostly monoterpenes, one phenylpropanoid (4-allyanisole), and one terpene derivative and its acetate ester (borneol and bornyl acetate) (hereafter monoterpenes) by comparing the mass spectra and retention times to those of known standards.

Monoterpenes were quantified using standard curves from four dilutions prepared from analytical standards of (+)-limonene (Chemical Purity: 99%), (-)-limonene (99%), α-terpinene (95%), (+)-α-pinene (98.5%), (-)-α-pinene (98%), myrcene (94%), 3-carene (98.5%), (+)-β-pinene (98%), (-)-β-pinene (99%), β-phellandrene (90%), (+)-camphene (90%), γ-terpinene (97%), terpinolene (90%), 4-allylanisole (98.5%), R-(+)-pulegone (97%), bornyl acetate (97%), (-)-borneol (99%), linalool (98%), ocimene (90%), and α-terpineol (90%). Concentrations were reported as “µg mg^-1 of dry weight (DW)” phloem tissue.

2.3.2. Diterpene Analysis

Diterpenes were extracted from 30 mg of lyophilized ground tissue in 1 mL of methanol. The mixture was vortexed for 30 s, then centrifuged for 10 min at 18,213 rcf at 4 °C. The supernatant was transferred to 2 mL GC vials and kept at -40 °C until analysis. Diterpene acid extractions were analyzed using Ultra-High Performance Liquid Chromatograph (UHPLC, 1290 Infinity, Agilent Tech.), fitted with an InfinityLab Poroshell 120 EC-C18 column (2.1 x 150 mm 1.9 µm, Agilent Tech.) and a Diode Array Detector (UV/Vis, 1290 DAD, Agilent Tech.). Gradient analyses were performed with a binary solvent system: 1.7% v/v acetic acid with ultra-pure water (A) and 100% methanol (HPLC grade) (B) and flowing at 0.3 mL min^-1. A 5 µL volume was used. The system started at 75% B for 1 min, then increased to 85% B over 9 min, held for 2 min, decreased to 75% B over 2 min and held at 3 min. The total run time was 17 min. Diterpene resin acids were quantified using wavelengths of 240, 268 and 282 nm as adapted from a multiwavelength detection approach from Kersten et al. [46]. Standard curves calculated from dilutions prepared of pimaric acid (80%), abietic acid (95%), dehydroabietic acid (99%), isopimaric acid (99%), levopimaric acid (95%), neoabietic acid (99%), palustric acid (92%) and sandarocopimaric acid (90%) were used to quantify diterpene acids. Pimaric and isopimaric acids were listed together due to coelution with abietic acid and cannot be quantified individually. Concentrations were reported as “µg mg^-1 of DW” phloem tissue.

2.3.3. Non-Structural Carbohydrate (Soluble Sugars and Starch) Analyses

All three soluble sugars (glucose, fructose and sucrose) were extracted from 30 mg of the lyophilized and ground phloem tissue. The sample was mixed with 1.3 mL of ultra-pure water in a 2 mL tube. A marble was placed in the opening of the tube. Tubes were placed above a boiling pot of water and enclosed in steam for 1 hr. A 0.5 mL aliquot was collected and centrifuged. Then, 0.4 mL of supernatant was placed in a new 2 mL tube with 1.0 mL HPLC-grade methanol. The extracts were then incubated at room temperature for one h. Then 0.5 mL of the extract was transferred to 2 mL GC vials and kept at -40 °C until analysis.

Glucose, fructose, and sucrose were analyzed with the UHPLC system that was fitted with an Infinity Lab Poroshell 120 HILIC-Z column (2.1 x 100 mm 2.7 µm, Agilent Tech.) as well as with an Evaporative Light Scattering Detector (ELSD, 1290 ELSD II, Agilent Tech.). A gradient analysis was performed using a binary solvent system of 0.034% v/v ammonium hydroxide to buffer ultra-pure water (A) as well as acetonitrile (HPLC-grade), also buffered with 0.034% v/v ammonium hydroxide (B). The flow rate was set to 0.2 mL min^-1. A 4 µL injection volume was used. The solvent system started at 90% B for 2 min, then decreased to 75% B over 3.5 min, where it was held for 1 min. It was then increased to 90% B over 2 min and held there for 0.5 min. The total run time was 9 min. The ELSD settings were set so that the nebulizer temperature was 25 °C, the evaporator tube temperature was 60°C, and the gas flow rate was 1.40 SLM. Dilution prepared from analytical standards of glucose (99%), fructose (99%) and sucrose (99.55%) were used to calculate standard curves to quantify soluble sugars.

For starch extraction, 25 mg of the lyophilized ground phloem tissue was used. Enzymatic digestions were used to convert starch into gluconate-6-phosphate, as described in Cale et al. [45]. As with the soluble sugar extraction, samples were held in steam for 1 hr. Samples were then vortexed for 30 sec. A 0.4 mL aliquot was transferred to a new 2 mL tube with 0.4 mL of alpha-amylase solution of 75 g enzyme (Sigma-Aldrich, St. Louis, MO, USA) in 100 mL ultra-pure water. This enzyme converts starch to large soluble sugars (e.g., maltose). The tubes were vortexed immediately for 30 sec and then incubated in a 50 °C bath for 16 hr. Following the incubation, tubes were removed from the water bath and then inverted twice. The extracts were then centrifuged at 18,213 rcf for 15 min, which removed undigested glucose-containing polymers (e.g., hemicellulose). After which 0.4 mL of the supernatant was transferred to a new 2 mL tube containing 0.4 amyloglucosidase solution consisting of 3 g enzyme (Sigma-Aldrich, St. Louis, MO, USA) in 60 mL sodium acetate buffer (0.1 M, pH 4.5). This enzyme converts maltose and other large sugars to glucose. The tubes were incubated in a 50 °C bath for 16 hr. Tubes were removed from the water bath and inverted twice. The extracts were centrifuged at 18,213 rcf for 15 min. Then, 0.4 mL of the extracts were transferred to 2 mL GC vials and stored at -40 °C until analysis.

We pipetted 0.020 mL of starch extract into a 96-well plate. Next, 0.2 mL of glucose hexokinase (Sigma Aldrich) - isomerase (phosphoglucose isomerase, Sigma Aldrich) was added to each well. This converted starch aliquots into gluconate-6-phosphate, as described in Cale et al. [45]. The well plate was shaken at room temperature for 45 min on an orbital shaker. Then, the amount of gluconate-6-phosphate was measured at an absorbance of 340 nm using the Synergy Microplate Reader H1 (BioTek, Winooski, VT, USA).

The starch concentration was quantified using two calibration curves: the first to estimate the glucose concentration from the amount of gluconate-6-phosphate measured and the second to estimate the starch concentration using the estimated glucose concentration. To do this, the sample concentration of total glucose was standardized by sample dry weight and extraction volume after the first calibration curve. Previously determined concentrations of soluble glucose and fructose (hexoses) were subtracted from the estimated total glucose concentration. A conversion factor was applied to account for potential quantification differences. The second calibration curve used dilutions of pure potato starch standards processed with the same enzymatic digestions as described above. The potential matrix effects were calculated to adjust the final concentrations accordingly.

2.4. Anatomical Defence and Growth Characteristics

Core and wedge samples were glued to wooden mounts and dried for two weeks at room temperature. Samples were then sanded with belt and hand sanders using progressively finer sandpaper. Samples were then scanned at 1,200 dpi to create high-resolution digital images. Annual growth (mm) was measured from pith to last year’s annual ring using ImageJ on all core samples [8]. Calendar years were assigned to each annual ring by visual cross-dating ring characteristics, such as width ratio and the colour contrast between earlywood and latewood per ring. This information was used for wedges to determine the year of death of trees killed by MPB. The earliest death recorded was considered the starting year of the outbreak in a particular site. Among the five sites, two sites experienced the outbreak in 2014, and the remaining three sites had an MPB outbreak in 2015. Therefore, the healthy trees were separated into two outbreak years: 2014 (n=24) and 2015 (n=36).

All samples for resin ducts were examined using ImageJ within a fixed sampling width of 10 mm [8]. To analyze resin ducts, resin duct production (number of resin ducts per 10 mm width per year on a ring [no. yr^-1]), resin duct size (mean size of resin ducts per 10 mm width per year on a ring [mm² yr^-1]), total resin duct area (sum of resin duct area per 10 mm width per year on a ring [mm² yr^-1]), resin duct density (total number of resin ducts per year divided by the ring area (10 mm * ring width) for a given year [no. mm² yr^-1], and relative resin duct area (percent area occupied by resin ducts per year within a ring area for a given year [% yr^-1]) was calculated for each 12 mm dia. core [13,45]. Ring width (mm yr^-1) and basal area increment (BAI) (mm² yr^-1) were used to represent tree radial growth. Tree radii and ring width data were used to calculate BAI, assuming that tree rings are concentric circles. Tree radii were calculated by dividing tree DBH by two.

2.5. Data Analysis

To study the long-term relationship between growth and defence mechanisms in healthy pine trees, we analyzed annual growth and resin duct characteristics across three distinct periods: two pre-outbreak (5-yr and 10-yr before the outbreak) and one post-outbreak (from outbreak onset until 2019 when the last annual ring was formed at the time of sampling) periods.

We conducted Pearson correlation tests to examine correlations between terpenes (21 monoterpenes, six diterpene resin acids, total monoterpenes, and total diterpenes), NSCs (three soluble sugars, starch, and total NSCs) and five anatomical defence traits and two radial growth characteristics. Analysis was performed using R (Version 4.0.5) [47] using Hmisc R package [48]. We selected BAI to represent tree radial growth as it strongly correlates with ring width (Pearson Correlation Coefficient > 0.96 in both the 2014- and 2015- outbreak years). In analyses, we considered overall as well as total and individual concentrations of constitutive monoterpenes and diterpenes.

2.5.1. Gradient Analysis

We employed an indirect gradient analysis to examine the impact of resin duct characteristics, BAI, tree age, and DBH on the overall monoterpenes, diterpenes, and NSCs before and after outbreaks. We utilized non-linear multidimensional scaling (NMDS) with Bray–Curtis dissimilarity to ordinate monoterpenes, diterpenes, and NSCs using the vegan R package [49]. NMDS results were saved as RData files to ensure reproducibility. Subsequently, we projected resin duct characteristics from pre- and post-outbreak periods onto NMDS axes as environmental gradients and visualized these relationships using ordination triplots created with the ggplot2 R package [50]. We followed the same steps when incorporating radial growth characteristics from pre- and post-outbreak periods and tree growth traits (age, DBH) into the NMDS axes for the 2014 and 2015 outbreak years.

2.5.2. Linear Mixed-Effect Models

To identify relationships among defence and growth characteristics, we employed linear mixed-effect models, treating the site as a random effect using the R package lme4 [51]. We visually assessed model assumptions, including normality and homogeneity of residuals, applying transformations and rescaling as needed. Resin duct density and the percent resin duct area were standardized per 1 mm² and accounted for BAI prior to calculation. Furthermore, we investigated the relationships between total monoterpenes or diterpenes and tree growth traits using linear mixed-effect models. When we observed a statistically significant relationship, we conducted further linear mixed-effect models with the site as a random effect to examine the specific relationship between individual monoterpenes or diterpenes and a particular predictor. In all cases, significance was determined at α = 0.05.

2.5.3. Comparative Analyses

We compared anatomical defences and radial growth responses across pre- and post-outbreak periods within the 2014 and 2015 outbreak years using linear mixed-effect models. Variability across sites was accounted for by nesting trees within sites. Furthermore, we examined the differences in tree age and DBH between the 2014 and 2015 outbreak years using linear mixed-effect models, where the site was considered a random effect to address any unexplained variances among trees.

3. Results

3.1. Do Historical Resin Duct Characteristics Relate to the Current Year’s Terpenes?

3.1.1. Monoterpenes

We observed several significant relationships between some of the resin duct characteristics and overall monoterpenes across both the 2014- and 2015-outbreak years (Figure 1, Table 1). In the 2014 outbreak, the relationship between the total monoterpenes and resin duct production and the total resin duct area was significant during all three periods. Resin duct size was also positively correlated with the total monoterpenes in the post-outbreak and 10-year post-outbreak period. In the 2015 outbreak, the total resin duct area and the resin duct size displayed significant positive relationships with total monoterpenes in all three time periods. Further testing of the relationship between individual monoterpenes and three resin duct predictors (resin duct production, total resin duct area, and individual resin duct size) yielded significant results. Individual monoterpene concentrations across sites are shown in Table S2. All three resin duct characteristics exhibited positive trends with several monoterpenes in the 2014 and 2015 outbreak years (Table 2). In the 2014 outbreak, we observed significant relationships between resin duct characteristics and 16 monoterpenes, while 15 monoterpenes displayed significant relationships with resin duct characteristics in the 2015 outbreak. Among the monoterpenes examined, several had significant relationships with at least one of the resin duct characteristics. Notably, (+)-limonene, (-)-α-pinene, (+)-β-pinene, β-phellandrene, myrcene, (+)-camphene, (-)-borneol and linalool correlated with most of the resin duct characteristics in both outbreak years (Table 2). Additionally, α-terpinene, (+)-α-pinene, 3-carene, γ-terpinene, terpinolene, (R)-(+)-pulegone, bornyl acetate, ocimene, and α-terpineol also exhibited significant relationships with at least one of the resin duct characteristics.

We did not observe any relationship between any of the resin duct characteristics and diterpene profiles or total diterpenes or individual diterpenes in either the pre- or post-outbreak periods within the 2014 and 2015 outbreak years, suggesting no association between phloem diterpenes and xylem anatomical defences (Figure 2).

Table 3. Results of linear mixed-effect models showing the relationship between total monoterpenes and total diterpenes, and the annual basal area increment (BAI) of lodgepole pine trees in Jasper National Park (Alberta, Canada) for different mountain pine beetle outbreak years and periods.

Outbreak Years	BAI	Total monoterpenes			Total diterpenes
		Coefficient	t-value	P-value	Coefficient	t-value	P-value
2014	Post-outbreak	532.4	2.7	0.015	-0.8	-0.5	0.608
	5-yr pre-outbreak	549.9	2.6	0.015	-0.8	-0.5	0.658
	10-yr pre-outbreak	665.1	3.2	0.004	-1.2	-0.7	0.516
2015	Post-outbreak	485.9	2.6	0.013	4.9	2.7	0.012
	5-yr pre-outbreak	686.8	3.7	0.002	4.9	2.4	0.024
	10-yr pre-outbreak	791.7	3.7	0.001	6.3	2.6	0.013

Significant P-values (α=0.05) were bolded and italicized.3.1.2. Diterpenes.

Table 3. Results of linear mixed-effect models showing the relationship between total monoterpenes and total diterpenes, and the annual basal area increment (BAI) of lodgepole pine trees in Jasper National Park (Alberta, Canada) for different mountain pine beetle outbreak years and periods.

Outbreak Years	BAI	Total monoterpenes			Total diterpenes
		Coefficient	t-value	P-value	Coefficient	t-value	P-value
2014	Post-outbreak	532.4	2.7	0.015	-0.8	-0.5	0.608
	5-yr pre-outbreak	549.9	2.6	0.015	-0.8	-0.5	0.658
	10-yr pre-outbreak	665.1	3.2	0.004	-1.2	-0.7	0.516
2015	Post-outbreak	485.9	2.6	0.013	4.9	2.7	0.012
	5-yr pre-outbreak	686.8	3.7	0.002	4.9	2.4	0.024
	10-yr pre-outbreak	791.7	3.7	0.001	6.3	2.6	0.013

Significant P-values (α=0.05) were bolded and italicized.3.1.2. Diterpenes.

3.2. Does NSCs Relate to the Terpenes?

In the 2014 and 2015 outbreaks, there were no significant relationships between any of the NSCs and overall monoterpenes (Figure S1, Table S3) or diterpenes (Figure S2, Table S3). However, when we tested the relationship between the total NSCs and the total monoterpenes or total diterpenes within the 2015 outbreak, the NSCs were significantly and positively related to the total diterpenes (coefficient = 0.701, P = 0.019) and the total monoterpenes (coefficient = 59.25, P = 0.048). We did not observe any significant relationships in the 2014 outbreak year.

3.3. Do the Radial Growth Characteristics Relate to the Terpenes?

Due to the strong correlation between BAI and the annual radial growth rate (ring width), we opted to use BAI as a representation of radial tree growth to account for tree size. Overall, monoterpenes are strongly associated with BAI across all pre- and post-outbreak periods (2014-outbreaK: r² values ranged from 0.414 to 0.516, P-values from 0.01 to 0.001; 2015-outbreak: r² values from 0.194 to 0.202, P-values from 0.047 to 0.033; Figure 3). Conversely, the relationship between diterpenes and BAI was marginally significant during the 10-yr pre-outbreak period in the 2015 outbreak (P = 0.05, Figure 4, Table S4).

Our analysis revealed strong positive relationships between individual monoterpenes and BAI (Table 4). BAI had a significantly positive relationship with 65% of individual monoterpenes in all pre-and post-outbreak periods in the 2014 outbreak and 50% of individual monoterpenes in the same periods in the 2015 outbreak. Specifically, (+)-limonene, (-)-α-pinene, myrcene, (+)-β-pinene, β-phellandrene, (+)-camphene, and linalool exhibited consistent positive relationships with BAI across all time periods in both outbreak years. However, this relationship appears to be stronger in the 2014 outbreak. Additionally, within the 5-yr and 10-yr pre-outbreak periods, the relationship between individual monoterpenes and BAI was more pronounced.

The analysis revealed a strong positive relationship between individual diterpenes and BAI across all pre- and post-outbreak periods, specifically within the 2015 outbreak (Table 4). This relationship extended to various diterpenes, including dehydroabietic, levopimaric, palustric, neoabietic, and abietic acids. In contrast, no significant relationship was observed between individual diterpenes and BAI within the 2014 outbreak.

3.4. Do tree DBH and Age Relate to the Terpenes?

Even though the trees sampled had similar DBH in 2014 and 2015 (2014: 22.19±0.59 cm; 2015: 20.76±0.32 cm), there was a notable difference in the mean age of trees between the two outbreak years (2014: 64.27±3.01; 2015: 80.42±3.85, P = 0.005). While DBH did not correlate with monoterpenes in either outbreak year, tree age exhibited a significant relationship with monoterpenes in both the 2014 and 2015 outbreaks (P₂₀₁₄ = 0.001, P_{2015 =} 0.017, Figure 5, Table S5). Tree age negatively correlated with the total (P₂₀₁₄ = 0.043, t-value: -2.2; P₂₀₁₅ = 0.009, t-value: -2.8; Figure 5) and individual (Table 5, Figure 5) monoterpenes. Tree age was linked to 45% of individual monoterpenes in the 2014 outbreak and 40% in the 2015 outbreak. Briefly, (+)-limonene, (-)-α-pinene, γ-terpinene, and linalool strongly correlated with tree age in both outbreak years.

The analysis indicated that tree DBH and age had no significant effects on diterpenes. Additionally, individual diterpenes did not exhibit any relationships with DBH or age (Figure S3) within the 2014 or 2015 outbreak years. Given these findings, further tests examining the relationships between individual diterpenes and tree age or DBH were not pursued, as no significant associations were detected.

3.5. Do NSCs Relate to Anatomical Defence Traits, Radial Growth Rate, and Tree Growth Traits?

The results from the indirect gradient analysis indicated that neither anatomical defences (Figure S4, Table S6) nor radial increment growth rate (BAI, Figure S5, Table S4) within the xylem had any substantial relationships with NSCs in the phloem for both the 2014 and 2015 outbreak years. Similarly, non-significant results were obtained when investigating the relationship between the total NSCs and each resin duct characteristic (Table 1) and BAI (Table 3) for both the 2014 and 2015 outbreaks. Additionally, neither tree age nor DBH appeared to impact the NSC profile, individual NSCs, or total NSCs (Figure S6).

3.6. Is There Any Trade-Off between Monoterpenes and Diterpenes?

We found no significant relationships between monoterpenes and diterpene in the 2014 or 2015 outbreaks (Figure S7, Table S7). Additionally, the relationship between the total monoterpenes and total diterpenes was insignificant in this study.

3.7. Additional Results

The correlation between anatomical defence and growth characteristics was significant across all time periods. The additional comparisons of anatomical defences and BAI among the post-outbreak, 5-yr pre-outbreak, and 10-yr pre-outbreak periods yielded significant results within both the 2014 (Figure 6) and 2015 (Figure 7) outbreaks. Overall, in both outbreaks, the resin duct characteristics exhibited greater size and relative area during the post-outbreak period compared to the 5-yr and 10-yr pre-outbreak periods. For instance, the annual resin duct size following the 2014 or 2015 outbreaks was larger by approximately 10-20% compared to both pre-outbreak periods. In the 2014 outbreak, the relative resin duct area after the outbreak was about 20% larger than in both the 5-yr and 10-yr pre-outbreak periods, even though BAI did not exhibit a significant response to the outbreak. In the 2015 outbreak, which comprised older pine trees, BAI during the 5-yr pre-outbreak period was 12% larger than in the 10-yr pre-outbreak period. However, no significant growth release was observed following the 2015 outbreak. Moreover, the total resin duct area during the post-outbreak period was 12% larger than during the 5-yr pre-outbreak period but similar to that observed in the 10-yr pre-outbreak period. Notably, resin duct production, density, and size displayed stronger responses following the 2014 and 2015 outbreaks than before. It's important to mention that including the site as a random effect in the linear mixed-effect models represented relatively high variance, which might contribute to the variations observed in terpenes and NSCs.

4. Discussions

We conducted a retrospective analysis of the constitutive monoterpenes and diterpenes, NSCs, and their interactions with anatomical defences and growth characteristics in surviving lodgepole pine trees in post-MPB-outbreak stands. We observed a strong correlation between post-outbreak concentrations of monoterpenes in the phloem and the historical resin duct characteristics in the xylem, along with a positive association between NSCs and total monoterpenes and diterpenes. Despite finding no evidence of growth release (BAI) in surviving pine trees, we observed strong positive relationships between BAI and several individual monoterpenes and diterpenes. Notably, tree age had a negative relationship with monoterpenes. Furthermore, the outbreak enhanced the resin duct characteristics in lodgepole pine. Our analyses, however, revealed no relationships between anatomical defence traits, BAI, tree age, size, and NSCs, nor did we find evidence supporting a trade-off between monoterpenes and diterpenes. These findings offer valuable insights into our broader understanding of conifer defences, particularly concerning how anatomical defences interact in response to outbreaking bark beetle species.

Our study showed strong correlations between monoterpenes and resin duct characteristics. Notably, overall monoterpenes were highly correlated with at least one of the five resin duct characteristics. For instance, the total monoterpenes were linked to the resin duct production, the resin duct size, and the total resin duct area. Furthermore, those monoterpenes that are crucial during MPB host colonization, such as (−)-α-pinene, β-phellandrene, (+)-limonene, and myrcene, had positive relationships with all resin duct characteristics, agreeing with our prior research findings [7,23,27,52,53]. Likewise, Kichas et al. [15] reported strong correlations between resin duct area, resin duct production and monoterpene concentrations in whitebark pine (P. albicaulis). In contrast, we found fewer relationships between diterpenes and resin duct characteristics, suggesting a weaker association between them. Similarly, Mason et al. [24] reported that anatomical defence traits of lodgepole pine differentially influenced the different classes of terpenes. Collectively, these findings highlight that the relationship between chemical and anatomical defences may depend on the classes of defence compounds [15]. This is in line with the conclusion from a meta-analysis, suggesting that trade-offs between monoterpenes and resin ducts are relatively rare across a wide range of plant species despite occasional occurrences in some systems [54]. Nonetheless, these results highlight the complementary role of anatomical and chemical defence structures in post-outbreak lodgepole pine stands, likely in response to the MPB outbreaks [7,8,9,22,24,40,42,43,55].

The second contribution of our study lies in the evidence regarding the positive relationship between total NSCs and total monoterpenes or total diterpenes within the 2015-outbreak. Although such relationships were not consistent among sites, this result suggests that the production of monoterpenes or diterpenes is linked to NSCs, aligning with previous research findings [32,33,34,35,56]. For instance, Mullin et al. [57] reported a positive relationship between starch concentration and total monoterpenes and diterpenes in lodgepole pine trees across varying elevations. Generally, it is well-documented that terpenoid defences require substantial energy demands [58,59,60] and deplete carbohydrate reserves [32,34,60]. However, the absence of a similar relationship between anatomical defences and NSCs suggests that our current understanding of the connections between sinks (defences) and sources (NSCs) in conifers is limited and necessitates further field testing [34,38].

One of the strongest relationships we observed is that both monoterpenes and diterpenes were positively correlated to tree radial growth, as measured by BAI. In fact, BAI exhibited positive relationships with up to 65% of individual monoterpenes, including (+)-limonene, (−)-α-pinene, myrcene, and β-phellandrene. Similarly, we observed a positive relationship between individual diterpenes, including dehydroabietic, levopimaric, palustric, neoabietic, and abietic acids, and BAI in the 2015 outbreak year. It is worth noting that the varying diterpene-BAI relationships between outbreaks in 2014 and 2015 could be attributed to differences in pine density, as the pine density in the earlier was nearly double that of the latter.

We observed negative correlations between tree age and monoterpenes in the 2014- and 2015-outbreak years. Specifically, the age of trees was linked to approximately 45% of individual monoterpenes in the 2014-outbreak and 40% in the 2015-outbreak. Notably, (+)-limonene, (−)-α-pinene, γ-terpinene, and linalool exhibited the strongest negative associations with tree age. These findings suggest that irrespective of tree size, tree age may be an important factor in determining the allocation of carbon resources to defence mechanisms. It appears that as trees mature, the allocation of resources to defences becomes less of a priority, whereas younger trees may invest more heavily in defence mechanisms. This could be attributed to the fact that younger trees are potentially more vulnerable to herbivory and thus prioritize defence allocation more than older trees, as supported by Swihart and Bryant [62] and Erbilgin and Colgan [63]. In the case of lodgepole pine, it was found that young trees in the range of 30–50 years old exhibited greater resistance to ophiostomoid fungi carried by MPB compared to both younger and older trees [64]. Overall, the relationship between the ontogeny of investment in defence versus growth in long-lived trees remains a complex and unresolved topic within the field of research [2,65].

The MPB outbreak influenced various physiological and growth traits of surviving lodgepole pine trees. Firstly, post-outbreak conditions improved resin duct characteristics compared to pre-outbreak conditions, with resin duct size, relative resin duct area, and total resin duct area all increasing. This suggests that increased investment in anatomical defences following MPB outbreaks, agreeing with the results of earlier studies. For instance, Zhao and Erbilgin [8] reported increased resin duct investment in surviving lodgepole pine trees. Resins serve as the initial line of defence for trees against bark beetle attacks by providing physical and chemical protection [20,30,66,67]. The number and size of xylem resin ducts in pine trees can predict resin production in response to wounding or bark beetle attacks [37,38,40]. Pines that survive MPB attacks tend to have more and larger xylem resin ducts [8,9,42,43,55] and increased pre-attack resin production [68]. Based on the anatomical defence and terpene duct characteristics observed in this and other studies, it would be expected that beetles would show a preference for larger pine trees.

Finally, our study did not provide evidence of growth release (BAI) in the surviving lodgepole pine trees. Furthermore, we found no connections between anatomical defence traits and BAI or evidence suggesting any notable relationships between monoterpenes and diterpenes. These findings contradict those of Kichas et al. [15], who reported significant correlations between BAI and resin duct characteristics in whitebark pine and lodgepole pine. We suspect that the enhanced anatomical defences observed after outbreaks may have required additional resources to support their development and maintenance rather than growth. Alternatively, other conifer species, such as spruce, and understory woody plant species like aspen, may have taken advantage of the available resources following the MPB outbreaks in this study, affecting the growth of surviving lodgepole pine trees. This suggests a complex interplay of factors influencing post-outbreak stand dynamics and resource allocation in these forest ecosystems, highlighting the need for further research to understand these intricate relationships fully.