Nitrogen isotope enrichment predicts growth response of Pinus radiata in New Zealand to nitrogen fertiliser addition

The fertiliser growth response of planted forests can vary due to differences in site-specific factors like climate and soil fertility. We identified when forest stands responded to a standard, single application of nitrogen (N) fertiliser and employed a machine learning random forest model to test the use of natural abundance stable isotopic N (δ15N) to predict site response. Pinus radiata growth response was calculated as the change in periodic annual increment of basal area (PAI BA) from replicated control and treatment (~ 200 kg N ha−1) plots within trials across New Zealand. Variables in the analysis were climate, silviculture, soil, and foliage chemical properties, including natural abundance δ15N values as integrators of historical patterns in N cycling. Our Random Forest model explained 78% of the variation in growth with tree age and the δ15N enrichment factor (δ15Nfoliage − δ15Nsoil) showing more than 50% relative importance to the model. Tree growth rates generally decreased with more negative δ15N enrichment factors. Growth response to N fertiliser was highly variable. If a response was going to occur, it was most likely within 1–3 years after fertiliser addition. The Random Forest model predicts that younger stands (< 15 years old) with the freedom to grow and sites with more negative δ15N isotopic enrichment factors will exhibit the biggest growth response to N fertiliser. Supporting the challenge of forest nutrient management, these findings provide a novel decision-support tool to guide the intensification of nutrient additions.


Introduction
As the global demand for wood and fibre increases, production forests have a key role to play in the growing bioeconomy (Marchetti et al. 2014;Payn et al. 2015). Increasing production to meet that demand can be achieved through a number of forest management options, including the use of mineral fertiliser to improve site nutrition and forest growth (Clinton 2018;Powers 1999). The addition of mineral fertiliser, in particular N, which is often limiting for forest growth (Davis et al. 2015;Fox et al. 2007;Littke et al. 2014), has been used to increase forest productivity at an operational scale (Chappell et al. 1991;Fox et al. 2007). However, the response of a forest stand to N fertiliser addition can be highly variable, resulting in low confidence in N fertiliser as a forest management option to increase production (Smaill and Clinton 2016;Sucre et al. 2008). The consequences of such a challenge are both economic inefficiency and potential ecological harm by facilitating N leaching. Thus, there is a growing call for improved precision silvicultural techniques, particularly around nutrient management (Rubilar et al. 2018).
Planted forests are an important mechanism for New Zealand to achieve their climate goals and desire to transition to a bioeconomy (Clinton 2018;NZFOA 2019). These forests are dominated by Pinus radiata D. Don (90% planted forest land area) and have a range in site fertility from low fertility sand dunes through to very fertile intensively developed pastoral soils (Beets et al. 2019;Garrett et al. 2022;Ross et al. 2009;Watt et al. 2008). The New Zealand forestry sector has traditionally managed forest nutrition exclusively to reduce nutrient deficiency (Davis et al. 2015;Mead and Gadgil 1978); however, there is now a greater interest in managing soil nutrition to achieve productivity gains (Clinton 2018;NZFOA 2019). The New Zealand forest industry currently 1 3 predicts a productivity response with N fertiliser addition using a basic grid of foliage N concentration deficiency and stand thinning or openness (Hunter 1982). An approach with soil variables is available, however is limited to expecting a productivity response at only 4 years post fertiliser addition and includes only a small sub-set of variables (Hunter et al. 1986). Both approaches presented (Hunter 1982;Hunter et al. 1986) are basic and limited by the data available at the time. Novel approaches are emerging that can transform decision-support tools, such as machine learning algorithms that can deal with the multiple variables and their complex nonlinear relationships (Morellos et al. 2016). Ability to predict a productivity response with the addition of N fertiliser at any site in any given year post fertiliser application is an important priority for the New Zealand forest industry to meet their aspirations to improve forest productivity sustainably (Smaill and Clinton 2016).
Though many studies have attempted to predict the response of forest stand productivity after N fertiliser (Hart et al. 1986;Hunter 1982;Hunter et al. 1986;Lea and Ballard 1982;Lim et al. 2015;Littke et al. 2014Littke et al. , 2017McNeil et al. 1988;Miller et al. 1989;Sucre et al. 2008;Turner et al. 1977), they are continually challenged by a lack of success, results that are too site-specific as to be extensible, or too demanding of complex parameterisation data to be useful. Thus, the "holy grail" in nutrient management decision support is something that is easy to measure, broadly transferable, and strongly predictive. We believe that natural abundance N isotopic signatures (δ 15 N) hold just such promise. Natural abundance δ 15 N values in soil and foliage are long-term integrators of ecosystem N cycling processes (Robinson 2001;Amundson et al. 2003;Craine et al. 2009). For example, plant δ 15 N has successfully been used to characterise forest ecosystem N cycle responses to global change drivers (BassiriRad et al. 2003;McLauchlan et al. 2007). Moreover, use of plant δ 15 N and the δ 15 N enrichment factor (EF), which is relative enrichment of foliage to soil (δ 15 N EF = δ 15 N Foliage − δ 15 N Soil ), have been shown to positively correlate with soil N transformations like net N mineralisation and net nitrification, connecting such patterns to specific N cycling processes with direct relevance to plant fertility and ecosystem N loss (Garten Jr and van Miegroet 1994). Additionally, the δ 15 N EF gives an insight into site-level N dynamics by integrating N cycling processes between soil and plant (Amundson et al. 2003;Cheng et al. 2010;Craine et al. 2009;Emmett et al. 1998;Fang et al. 2011;Garten Jr and van Miegroet 1994). Though the application of this characteristic to predicting plant or ecosystem response to N fertiliser addition has never been attempted, the observations referenced above suggest that δ 15 N profiling should be precisely the type of integrative ecosystem characteristic that predicts how a forest will deal with added N. In fact, preliminary results have demonstrated that the δ 15 N EF is useful at predicting productivity response to N fertiliser addition across multiple species (e.g., Douglas-fir and loblolly pine; (Lorentz 2013)). Thus, our goal here is explore the potential of site-level measures of δ 15 N signatures to similarly predict the response of P. radiata, the predominant planted forest species in New Zealand, a system in which the technique has never been employed, but where such approaches at assessing plant-soil fertility relationships are increasingly called for (Smaill and Clinton 2016).
The objective of this study was to explore predicting a growth response in planted P. radiata forest stands, in New Zealand, after the addition of N fertiliser. We will explicitly test the hypothesis that sites with more negative δ 15 N EF signatures will have a larger tree growth response to N fertiliser. To test this hypothesis, we used non-parametric Random Forest (RF) machine learning modelling to evaluate the relative importance of foliage and soil δ 15 N values against other, more commonly measured ecosystem characteristics.

Study sites
Historical and current N fertiliser trials in New Zealand were used in this study. Criteria required for the inclusion in the study were: tree species was P. radiata, replicated trial design with a control and N fertiliser treatment, tree age over 3 years old when fertiliser applied, N fertiliser addition of a known amount and year of application, tree measurements of diameter at breast height (DBH; 1.4-m stem height) and tree height post fertiliser application, an archived soil sample (0-10 cm) and foliage sample (1-year-old foliage) collected before fertiliser application or from the control treatment. These criteria were used to ensure enough data from each site to properly reflect the site conditions. A total of 18 sites were identified (Fig. 1). Twelve sites from Scion's (New Zealand Forest Research Institute Ltd) historical and current experimental trials fulfilled all criteria: 5 sites within the FR467 trial series, 5 sites within the FR561 trial series, and 2 sites RO1843 and WN379 (Fig. 1). An additional 6 sites were identified that matched the criteria except for an archived soil sample (trial ID = AK976 (2 sites), RO1889, RO1083 (2 sites), and RO1818).
All 18 sites were used to test when a forest stand will respond to a single application of N fertiliser (150 -207 kg N ha −1 ). Eight of the trial sites (AK976 (2 sites), FR467 (5 sites), and RO1889) also had plots with an N + P (mixed range of 80 -400 kg N ha −1 and 40 -200 kg P ha −1 ) and P (75 kg P ha −1 ) single fertiliser additions, and these two treatments were also included in the analysis of when a forest stand will respond to a single application of fertiliser. Only the 12 sites that matched all criteria were used in the analysis of predicting a growth response after a single application of N fertiliser (200 -207 kg N ha −1 ). The range in site soil and foliage properties, plus tree age at the time of N fertiliser addition for the data sub-set used to predict a growth response, are listed in Table 1 (additional detail in Supplementary Table S1).

Sample preparation and analysis
Foliage samples were previously oven dried at 70 °C to constant weight and ground to < 1 mm prior to chemical analysis. Mineral soil samples were previously air-dried (< 40 °C) and then sieved to retain the < 2-mm fraction for Fig. 1 Distribution of trials included in the analysis of 1) testing when a forest stand will respond to a single application of fertiliser (using all trials shown), and 2) predicting a growth response ( 1 3 chemical analysis. Results from past chemical analysis were used, which included total carbon (C) and total N (measured by LECO Trumac CN (modified Dumas)), soil acidity pH (measured by electronic probe, 1:2.5 dH2O), foliage total phosphorus (P) using wet digestion (Nicholson 1984) for historic sites (RO1843 and WN379), and foliage total P measured by ICP-MS at all other sites. All chemical analyses are reported at 104 °C oven dry basis. New chemical analysis was undertaken for any samples without total C and total N, and all samples were analysed for δ 15 N natural abundance using a continuous flow isotope ratio mass spectrometer (IsoPrime 100 EA-IRMS, Isoprime © Ltd., Manchester, UK). Soil C/N ratio (total C / total N) and site δ 15 N enrichment factor (δ 15 N foliage − δ 15 N soil ) were calculated.

Data and statistical analysis
Published studies that have modelled a productivity response with fertiliser addition in production forests have used tree growth metrics like basal area, volume, periodic annual increment, and/or site index, and fit these to a regression model (Carter et al. 1998;Hunter 1982;Hunter et al. 1986;Littke et al. 2014Littke et al. , 2017McNeil et al. 1988;Sucre et al. 2008); however, they are ultimately constrained because the number of predictor variable can vastly outnumber observations. Machine learning algorithms, such as the widely used RF modelling (Breiman 2001), can overcome this limitation by enabling a larger number of covariates to be used as they are able to deal with complex nonlinear relationships between the predictor and response variables (Morellos et al. 2016). There have been increasing applications of machine learning approaches in forest research, including tree species and patch distributions (Iverson et al. 2008;Veldhuis et al. 2017), tree species richness and C storage (Lautenbach et al. 2017), forest allometric scaling relationships (Duncanson et al. 2015), and wind damage (Moore and Lin 2019). Thus, following other contemporary studies on forest responses to fertiliser addition (Littke et al. 2014(Littke et al. , 2017, we will employ these techniques as discussed below.

Predicting a growth response
Predicting a growth response after N fertiliser addition was investigated using a sub-set of the data that included all criteria (12 sites). Only the N fertiliser addition treatment was selected as this treatment had a good level of tree measurements post fertiliser addition for all 12 sites compared to the N + P treatment which had limited data (only FR467 was available for N + P analysis). Tree growth was represented by the change in periodic annual increment of basal area (PAI BA m 2 ha −1 year −1 ) at the plot level. To quantitatively evaluate the growth response and the importance of different features (site, silvicultural, environmental, and geochemical) after the addition of fertiliser, we employed a non-parametric machine learning approach called Random Forest (RF). This approach was chosen because of the complexity of our data (different experimental designs, low number of replicates, asymmetrical experimental design, and variable environmental factors), which made traditional statistical models inadequate. RF is an ensemble method of machine learning, which is based on averaging the prediction or taking the majority vote of a large group of classification and regression trees (CART) learners (Breiman 2001). Each CART is a tree prediction model fitted using random samples and variables, drawn from same distribution, from the original dataset. Detailed methods behind the Random Forest model are held in the Supplementary.
Features used in our RF model as predictor variables were from four groups: site, silvicultural, climate, and geochemical features (Table 2). Climate feature variables were sourced from the climate data of the nearest climate station in NIWA's Virtual Climate station Network (VCSN) (Tait and Turner 2005) to each plot with its location and time period. The VCSN data are estimates of several climate variables interpolated on a regular 5-km grid covering the whole of New Zealand from 1972 to present. The original data was recorded in a daily form, which did not fit our modelling purpose, so we restructured the data into monthly, seasonal, and annual forms. The annual climate dataset was used for modelling due to a balance between computational economy and model performance. There were 24 climate-related predictor variables. For those plots 1 3 planted and measured before 1972, we used 5-year climate data from 1972 to 1976 as a proxy.
To evaluate and select candidate predictor variables used in the RF model, we employed the Boruta algorithm (Kursa and Rudnicki 2010). The Boruta algorithm is a wrapper built around the RF algorithm, which provides a statistically grounded way for automatically selecting variables on a given dataset. According to Boruta evaluations, all the forty-eight candidate variables were significant for model improvement based on their importance ( Supplementary Fig. S1) and were used in our RF model as predictor variables. All analyses and machine learning models were implemented in the R 3.5.3 platform (R Core Team 2017). The RF models were implemented, trained, and validated by using R package randomForest (Liaw and Wiener 2002) and caret (Kuhn 2008). Variables most associated with tree growth PAI BA were determined by > 50% relative importance in the RF model (Duncanson et al. 2015). Partial dependence plots for each predictor variable were also generated to support interpretation of the model predictions.

Growth response
Duration of tree growth response to N fertiliser was investigated using the full dataset (18 trial sites) and all treatments: control treatment and N addition at all sites, and N + P and P addition at 8 sites. Tree growth was represented by the change in periodic annual increment of basal area (PAI BA m 2 ha −1 year −1 ), a representation of the annual growth rate for the year in question, at the plot level, and was determined using measured plot area and tree stocking. (1) There was a large variation among trial sites and plots (i.e., climate, soil, experimental design, replication, and silvicultural management). Therefore, to test for a significant response, we used a relative response index (RRI) to measure the growth response (PAI BA) of each plot at different years after fertiliser addition.
where X 0 is the mean performance (PAI BA) of control plots (no fertiliser addition) within a site and X 1 is the performance of treated plots (with fertiliser addition). RRI is symmetric and ranges from − 1 to 1. RRI > 0, positive response; RRI = 0, no response; RRI < 0, negative response. For each treatment, one tailed t tests were used to determine if plots' response (RRIs) is larger than 0 (i.e., only positive responses were counted). In the analysis, trials and treatments are considered independent experiments. We used Fisher's combined probability test on p values of RRI for each year post fertiliser addition.

Results
The twelve trial sites included in our RF model covered a large range in the on-site variables measured including δ 15 N EF which ranged from − 0.8 to − 8.4 (Table 1). The tree age at the time of fertiliser addition averaged 15 years with a range between 9 and 22 years of age (Table 1).
Overall, the best trained RF model was able to explain 78% of the variation in the productivity using PAI BA (R 2 = 0.778), showing a reasonable prediction performance according to the repeated cross-validation. Based on the RMSE value, this best trained RF model performs consistently well and has the smallest RMSE (0.479). In the RF T-mean (mean annual temperature), T-min (minimal annual temperature), T-max (maximal annual temperature), Rain-mean (mean annual precipitation), Rain-min (minimal annual precipitation), Rain-max (maximal annual precipitation), Rad-mean (mean annual solar radiation), Rad-min (minimal annual solar radiation), Rad-max (maximal annual solar radiation), PET-mean (mean annual potential evapotranspiration), PET-min (minimal annual potential evapotranspiration), PET-max (maximal annual potential evapotranspiration), RH-mean (mean annual relative humidity), RH-min (minimal annual relative humidity), RH-max (maximal annual relative humidity), SoilM-mean (estimated mean annual soil moisture), SoilM-min (estimated minimal annual soil moisture), SoilM-max (estimated maximal annual soil moisture), VP-mean (mean annual water vapour pressure), VP-min (minimal annual water vapour pressure), VP-max (maximal annual water vapour pressure), Frost-mean (mean annual frost days), Frost-min (minimal annual frost days), Frost-max (maximal annual frost days) model, tree age at the time of measurement (AGE-YEAR) and the enrichment factor (EF = δ 15 N Foliage − δ 15 N Soil ) were the two most important variables (> 50% relative importance) associated with tree growth PAI BA (Fig. 2). Partial dependence plots were used to illustrate marginal effect of the important predictor variables and the temporal changes of tree growth in response to the control and N fertiliser treatment ( Fig. 3 and single plots Supplementary Fig. S2). According to the RF model, P. radiata PAI BA was greater in younger stands < 15 years of age ( Fig. 3a) and greater where the δ 15 N EF was less negative and increasing sharply in a step change when δ 15 N EF was greater than − 2‰ (Fig. 3b). Partial plots showed PAI BA increased with N fertiliser treatment compared to the control (Fig. 3, and Supplementary Fig. S2). A response in PAI BA after fertiliser application was observed in trees younger than 20 years (Fig. 3a). Application of N fertiliser increased the PAI BA across the range of δ 15 N EF values (Figs. 3b and 4a). A fertiliser application of 200 kg N ha −1 showed greater growth response when the δ 15 N EF was more negative (Fig. 4b).
Results from the RF model showed greater growth responses within the first 3 years after N fertiliser addition (Fig. 3c). This finding is supported by the significant (P < 0.05) PAI BA growth responses in P. radiata that were observed within the first 3 years after N fertiliser addition using the full dataset (18 sites) (Fig. 5). Tree growth response to the N fertiliser addition for these 18 sites varied largely, with both positive and negative relative responses in PAI BA (Fig. 5). The addition of N + P fertiliser resulted in a longer period of growth response out to 4 years and 8 years after fertiliser addition (Fig. 5), while the addition of only P fertiliser resulted in no significant (P > 0.05) growth response (data not shown).

When will Pinus radiata respond to N fertiliser addition?
One of the challenges in managing the fertility of a production forest is evaluating the time course and duration of response to fertiliser addition. Previous research in New Zealand on predicting the growth response of P. radiata after N fertiliser addition of 200 kg N ha −1 had been limited by only reporting at 4 years post fertiliser addition (Hunter et al. 1986). The time period of 4 years selected by Hunter et al. (1986) was informed by the growth response from only N deficient sites (e.g., sand dune forests) where a response to N fertiliser was previously reported (Hunter 1982;Hunter and Hoy 1983). Our study demonstrated that fertiliser application rates of approximately 200 kg N ha −1 are more likely to result in a positive growth response in the first 3 years following application (Figs. 4c and 5). However, our results also showed that the addition of N + P fertiliser to a site may extend the time frame for a positive response (Fig. 5), highlighting the importance of considering other possible co-limiting nutrients, rather than focusing exclusively on N (Carter et al. 1998;Hunter et al. 1986;Sucre et al. 2008;Vadeboncoeur 2010).
In addition to the duration of the N fertiliser growth response, the RF model demonstrated that tree age was the most important variable for predicting a productivity Fig. 2 Relative importance of predictor variables in RF model for PAI BA. The importance of predictor variables was normalised by the maximum variable importance for the RF model to produce standardised measures of variable importance. YAF stands for year after fertiliser addition 1 3 response with N fertiliser addition. Our RF model showed that trees less than 15 years of age (mid-rotation age; 29-year average rotation length (MPI 2019)) have the greatest PAI BA regardless of treatment (Fig. 3a). Experimental trials included in our RF study were between 9 and 22 years of age when fertiliser was applied (Table 1), and a marginally greater growth rate in N-application stands was visible up to 20 years of age (Fig. 3a). Therefore, as there was less growth response of trees older than 15 years of age to N fertiliser addition, the addition of fertiliser may not be economically efficient. Studies have identified that the best response from fertiliser addition is when the trees have "free-to-grow" conditions, which are created through stand thinning (Carter et al. 1998;Sucre et al. 2008). For example, Hunter et al. (1986) observed that large basal area responses to N fertiliser addition in New Zealand P. radiata occurred where stands were less than 10 years old and soils were N-poor.
However, overall there are inconsistent results from research on P. radiata stand growth response to N fertiliser in New Zealand (Hunter et al. 1986;West 1998;Woollons and Will 1975). The dataset used for our RF model had a range of silvicultural treatments, including thinning and no thinning, and to test the effects of thinning and fertiliser application on growth response would require more sufficient datasets from rigorous experiments. However, our results confirm that there was a decline in PA BAI after age 15 years and support fertiliser application aimed at boosting tree growth during rapid-growth-rate "free-to-grow" periods.

Use of N stable isotopes predicting a productivity response with N fertiliser addition
In this study, δ 15 N EF was the second most important variable in the RF model for predicting a productivity  response with N fertiliser addition. Our hypothesis was shown to be correct, that a tree growth response to N fertiliser addition is more likely on sites with a more negative δ 15 N EF. The δ 15 N EF provides a single value of site N status that integrates historical climate and N cycling activity (Garten Jr and van Miegroet 1994) and combined with non-parametric RF machine learning modelling has proven to be a strong variable to include in New Zealand P. radiata nutrition analysis. Moreover, there is merit in the further exploration of the utility of the δ 15 N EF in forest nutrient management with thesis research by Lorentz (2013) showing that the δ 15 N EF was a useful metric in predicting Douglas-fir and loblolly pine productivity response to N fertiliser addition in the Pacific Northwest and southeastern USA, respectively. Other studies also found site fertility variables important in predicting a productivity response such as surface soil and forest floor C/N ratios (Littke et al. 2014), soil C andN concentrations (Sucre et al. 2008), soil total N and available phosphorus (Hunter et al. 1986), and soil mineralisable-N (Carter et al. 1998). Our results further emphasise the importance of considering site fertility, including some metric of N cycling (i.e., δ 15 N) to predict N fertiliser growth responses.
Soil and foliage enrichment factors in forests have been shown to be related to overall rates of soil N dynamics, cycling, and N loss (Cheng et al. 2010;Emmett et al. 1998;Fang et al. 2011;Garten Jr and van Miegroet 1994). It could be expected that sites with a more negative δ 15 N EF would have lower rates of net N mineralisation and net nitrification (Garten Jr and van Miegroet 1994) and therefore be more N limited and conservative with their N cycle. We found that along the gradient of δ 15 N EF, growth rates generally decreased with more negative enrichment factors, in both the control and fertiliser addition treatments (Fig. 3b). With the addition of 200 kg N ha −1 fertiliser, our results showed the growth response was greatest at sites with a more negative δ 15 N EF (Fig. 4b). Both of these observations are consistent with the conceptual model described by Gurmesa et al. (2022), in which N-limited sites would have the most negative enrichment factors.
Where sites have low rates of N mineralisation and net nitrification N can be supplied by mycorrhizal fungi, which may transfer strongly δ 15 N-depleted N to the trees. Such N-limited sites would then be expected to exhibit strong growth responses to added N (as seen in this study; Fig. 4). As N availability increases, and dependence on mycorrhizal fungi decreases, enrichment factors become less negative, as trees shift to direct uptake of soil mineral N (Gurmesa et al. 2022). Nitrogen inputs and exports from other sources are also important to consider in understanding a site's δ 15 N dynamics (Högberg 1997). Inputs of N from atmospheric deposition in New Zealand are relatively low with values ranging from 3 to 9 kg N ha −1 year −1 (Parfitt et al. 2006) and therefore unlikely to have a strong impact on the soil and foliage δ 15 N at our sites. Inputs of N from N-fixing weeds or understory species may also be possible, with on average 30 kg N ha −1 year −1 reported where N-fixing species are present in New Zealand planted forests (Parfitt et al. 2006). Fires which can reduce soil N were used in historical forest clearance in New Zealand, although, fires in planted forests have been low with only 2% of planted forest land area impacted (McIntosh et al. 2005;Pearce et al. 2008). We are uncertain of any past fire events or presence of N-fixing plants at our sites used in the RF model; therefore, N input from fixation or loss of N through fire, and their resulting impacts on soil and foliage δ 15 N are unknown (Högberg 1997). We observed a step change increase in productivity as the enrichment factor increased beyond − 2‰ (Fig. 4a); moreover, the growth response to fertiliser also increased in contrast to the general decrease in fertiliser response seen up until that point (Fig. 4b). It is possible that this step change in δ 15 N EF at − 2‰ may be a threshold indicating a different mechanistic shift in N acquisition strategy. For example, it could indicate a different degree of mycorrhizal association or it could be an integration of multiple factors including ecosystem-scale additions and losses. However, more data is required to explore this possible threshold and its mechanistic underpinnings. We also observed that for enrichment factors more negative than − 5‰ there appears to be no change in the gain in productivity with 200 kg N ha −1 fertiliser addition (Fig. 4b). This indicates that there was a cutoff in expected growth when sites have very negative enrichment factors, possibly because these very negative δ 15 N EF sites require more than a single addition of 200 kg N ha −1 fertiliser to meet tree N demand (e.g., Smaill et al. (2011)). This threshold may alternately be reflecting other nutrient limitations that are constraining any further growth response to the added N (Vadeboncoeur 2010). Overall, the δ 15 N EF moves us towards a quantifiable decision support tool for forest nutrient management decisions. However, beyond the thresholds described above the δ 15 N EF was sufficiently vague with regard to the mechanistic drivers behind the result to allow for the rationalisation of alternative nutrient management strategies.

Implication for New Zealand forestry
In New Zealand, planted forests are located on sites with a wide range in fertility, derived from differences in soil parent material and past land-use, such as pastoral farming history and fertility enhancement through fertiliser addition (Beets et al. 2019;Garrett et al. 2022;Ross et al. 2009;Watt et al. 2008). Currently, the New Zealand P. radiata forest industry relies on the use of foliage N concentration and target ranges to identify N deficient stands, which would benefit from N fertiliser application (Hunter 1982). However, this strategy results in highly variable growth responses (Smaill and Clinton 2016). The 4-year basal area response model by Hunter et al. (1986) was an advancement on using foliar N concentration and also includes soil nutrient variables (total soil N and available phosphorus), soil clay percent, tree age, and simple silviculture variables (pruning and thinning), which together explained 66% of the variation in growth response.
Our RF model is a novel approach towards modelling the complexity of forest ecosystems and explains 78% of the variation in forest stand growth variation as a result of N fertiliser addition. Our model predicts annual incremental growth metrics of the forest stand basal area (PAI BA) which allows for forest management on a more detailed scale to understand if a forest stand will respond more in year one compared to another year. Moreover, in developing the RF model we were able to include more complex data to inform the prediction, for example, detailed climate variables were included that were both site and time specific. Moreover, we have demonstrated that the pre-fertiliser δ 15 N EF was a strong variable in predicting if a site will respond to N fertiliser. Our RF model's performance was good (R 2 = 0.778), with a dataset from 12 sites, which is a relatively small training dataset for machine learning modelling. The RF model could be improved with more training data and should be validated with further independent datasets. Models that can predict a productivity response of a stand with pre-fertiliser site variables will enable focused fertiliser applications, targeted at achieving the greatest overall boost in productivity within a forest or a set of forests, to support forest productivity enhancement.
Our RF model predicts a productivity response in PAI BA, which can be used to support site-specific nutrient models in predicting a growth response with N fertiliser addition. An example is the New Zealand Nutrient Balance Model (NuBalM) (Smaill et al. 2011), which describes the annual N demand and supply within a forest stand. This model is currently unable to predict a growth response with N fertiliser addition. Linking growth response predictions of the RF model to NuBalM would enable increasingly accurate N balance and cycling predictions over multiple rotations. Furthermore, this approach will provide a basis for sitespecific recommendations of mineral N fertiliser necessary to improve New Zealand planted forest productivity.

Conclusion
Measurement of the soil and foliage δ 15 N and use in a sitespecific δ 15 N EF has been shown in our study to be a strong variable in understanding forest N dynamics in response to newly added N. Using a machine learning RF model to deal with complex site-specific data, the model showed tree age and the N isotopic EF to be the most important variables affecting P. radiata growth in New Zealand. Our RF model can be used to predict site-specific annual tree growth response with the addition of 200 kg ha −1 of N fertiliser. A growth response of P. radiata to 200 kg ha −1 of N would be expected to occur within the first 3 years after fertiliser addition; however, our study showed both positive and negative growth responses within those first 3 years can also occur. Tree growth rates can be expected to be lower where a site δ 15 N EF is more negative, and with the addition of N fertiliser, these sites, with more negative δ 15 N EF, are more likely to have a growth response. The ability to make site-specific predictions will enable forest managers to 1 3 selectively target sites for productivity improvement. Future development of the RF model through more observations could significantly improve the model's ability to predict a productivity response with N fertiliser addition across a larger range of climate and soils. Further studies should then include validation of the model with an independent dataset to ensure the robustness of predictions. There is also potential for the wider application of a site's δ 15 N EF in predicting a forest productivity response to N fertiliser addition in range of forestry systems.

Acknowledgements
We thank the forestry mangers who provided sites for the fertiliser trials and the Scion staff for sample sourcing and result discussions.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.