Quantitative Genetic Variation in Bark Stripping of Pinus radiata

: Bark stripping by mammals is a major problem for conifer forestry worldwide. In Australia, bark stripping in the exotic plantations of Pinus radiata is mainly caused by native marsupials. As a sustainable management option, we explored the extent to which natural variation in the susceptibility of P. radiata is under genetic control and is thus amenable to genetic improvement. Bark stripping was assessed at ages four and ﬁve years in two sister trials comprising 101 and 138 open-pollinated half-sib families. A third younger trial comprising 74 full-sib control-pollinated families was assessed at two and three years after planting. Signiﬁcant additive genetic variation in bark stripping was demonstrated in all trials, with narrow-sense heritability estimates between 0.06 and 0.14. Within sites, the amount of additive genetic variation detected increased with the level of bark stripping. When strongly expressed across the two sister trials, the genetic signal was stable (i.e., there was little genotype × environment interaction). No signiﬁcant non-additive e ﬀ ect (speciﬁc combining ability e ﬀ ect) on bark stripping was detected in the full-sib family trial, where it was estimated that up to 22.1% reduction in bark stripping might be achieved by selecting 20% of the less susceptible families. Physical traits that were genetically correlated, and likely inﬂuenced the amount of bark removed from the trees by the marsupials, appeared to depend upon tree age. In the older trials, these traits included bark features (presence of rough bark, rough bark height, and bark thickness), whereas in the younger trial where rough bark was not developed, it was the presence of obstructive branches or needles on the stem. In the younger trial, a positive genetic correlation between prior height and bark stripping was detected, suggesting that initially faster growing trees exhibit more bark stripping than slower growing trees but later develop rough bark faster and became less susceptible. While the presence of unexplained genetic variation after accounting for these physical factors suggests that other explanatory plant traits may be involved, such as chemical traits, overall the results indicate that selection for reduced susceptibility is possible, with potential genetic gains for deployment and breeding.


Introduction
Herbivores are important determinants of plant productivity in managed and natural plant systems [1]. In managed conifer forests, browsing may have significant economic impacts through between the ages of 1 and 6 years, has become the most important pest problem [6,38,46], affecting up to 40% of the Tasmanian plantations, with up to 80% of trees affected in some plantations [6]. Given that the financial costs of managing mammalian herbivores through fencing and culling in P. radiata plantations are high [47], the potential for genetic variation in resistance to bark stripping to sustainably reduce damage is of interest.
This study aimed to: (1) determine the extent to which bark stripping is under additive genetic control and if genetic differences are stable across sites and tree age; (2) determine the genetic correlation between the level of bark stripping, growth, stem and bark traits; (3) estimate possible genetic gains in reducing bark stripping damage from the field-based selection of the least damaged families.

Family Trials
The field trials used in this study were established in Tasmania ( Figure 1) by Timberlands Pacific Pty Ltd. with seeds from the Radiata Pine Breeding Company (New Zealand). Three trials were studied-Beulah, Payanna, and Wilmot ( Table 1). The progeny planted in Beulah and Payanna consisted of fourth-generation half-sib (open-pollinated; OP) families from selections in the breeding program, while that of Wilmot was established from third-generation full-sib (cross-pollinated; CP) families. The Wilmot CP families (n = 74 families) were derived from 55 parents and 54 unique grandparents. The OP families at Beulah (n = 101 families) and Payanna (n = 138 families) were assumed to be half-sibs [48] and were derived from 101 and 138 mothers and 194 and 195 grandparents, respectively, highlighting the high pedigree diversity of the populations in the study and that the results were not biased by a few founder ancestors ( Table 1). The progeny were from parents that were selected for vigor, stem form, wood properties, and branch characteristics [39,40].  Table 1. Details of the three Pinus radiata genetic trials used for the study indicating the location and date of planting as well as the trial design, pedigree information, and assessment information. All trials were replicated in a randomized incomplete block design, with open-pollinated (OP) or cross-pollinated (CP) families represented as single tree plots within each block. For each trial, the number of replicates and incomplete blocks as well as the families, parents, grandparents, and total number of trees are indicated. The total number of trees denotes the trees that were alive at the time of assessment. Selected traits (see results) were measured multiple times and the age at which the assessments were performed is shown.

(5 years)
Forests 2020, 11, 1356 5 of 24 The trials planted at Beulah and Payanna were planted in the same year, with 98 common families plus 3 families which were unique to Beulah and 40 families which were unique to Payanna. The genotype by environment interaction (G×E) was tested based on these two trials. The younger Wilmot trial did not have parents in common with those represented in the other two trials but had some common ancestry deeper in the pedigree. All the trials were replicated in a randomized incomplete block design, with families represented as single tree plots within each block. Seedlings were raised in pots and were planted in rows, at a spacing of 3 m by 3 m. Since the trials were established in clear-cut coupes formerly dominated by Pinus radiata, the remaining post-harvest debris was gathered to form wind rows which separated different sets of blocks within each trial.
In Tasmania, the major bark stripping pest of P. radiata-the Bennett's wallaby (Macropus rufogriseus) [46]-varies in population density depending on location. The amount of alternative wallaby food may also vary in the different sites, although the presence of other food may not reduce bark stripping [49]. The Beulah and Wilmot field trials were situated in mid-north Tasmania, where the density of the Bennett's wallaby is estimated at approximately 32.0 animals/km 2 [50]. The Payanna field trial was situated in the north-east, where the estimated density of the Bennett's wallaby is approximately 26.9 animals/km 2 [50]. Apart from operational culling [50], the field trials at Beulah and Payanna were freely accessible by the animals. The Wilmot site was, however, fenced for the first two years of growth. After two years, the gates to the trial were opened during winter for approximately 2 months to allow the marsupial herbivores access to a subset of 20 of 26 replicates (6 replicates spread throughout the trial remained independently fenced). The gates were closed to stop browsing when browsing was evident across all 20 replicates. However, after the first bark stripping assessment animals breached the fence and accessed the trees. Further bark stripping was also scored.

Assessment of Bark Stripping Damage and Related Traits
Bark stripping damage on the stems of the P. radiata trees was recorded visually on individual plants in the field trials. At Beulah and Payanna, the damage was scored on an ordered categorical scale, assigning zero (0) to plants with no evidence of bark stripping, 1 = < 20% of the circumference stripped, 2 = 20-50% of the circumference stripped, and 3 = 50-100% of the circumference stripped. At Wilmot, more categories were included-i.e., 3 = 50-75%, 4 = > 75%, 5 = 100% damage (completely ring barked). Except for the scores 0 and 100, the remaining scores were converted to class mid-point values for the final analyses. Other tree traits including tree height, basal diameter (BD), diameter at breast height (DBH = 1.3 m), bark thickness, the presence/absence of rough bark, stem access, and survival were also assessed (Tables 2-4). Rough bark was assessed by the visual inspection of the presence or absence of bark fissures (e.g., Figure 2), and, when present, rough bark height was measured as the height from the ground to the top of the rough bark on the side of the stem with the highest rough bark cover. Bark thickness was estimated at 1.3 m with a custom-made bark probe and was not corrected for DBH. Stem access describes the presence of stem needles and obstructive branches in the first 1 m of the stem, which may prevent ready access to the bark. This trait was subjectively assessed in a categorical decile scale, where 0 = all stem covered with needles and branches, 10 = up to 10 cm covered, 20 = up to 20 cm covered, until 100 where no needles or branches were found within 1 m. Stem access was important only in the younger Wilmot trial. In the older trials, the stem needles were dead and few trees had sparse branches on the stem within 1 m from the ground. Tree survival at all trials was assessed at the time of the first bark stripping assessment. Some variables were assessed multiple times (Tables 2-4). The first assessment at Payanna and Beulah was made at 4 years and then at 5 years after planting. At the Wilmot trial, the first, second, and third assessments were made at 2, 3, and 5 years after planting, respectively. For bark stripping, recurrent assessments included old damage, except where the tree had completely healed with no clear signs of earlier damage, and, therefore, the two scores are not independent.

Linear Models
All the univariate and bivariate linear mixed models were run in ASReml-R [51]. First, univariate models were fitted to generate variance components that were used to test the presence of additive genetic variation. The general linear mixed model is represented as: where y is the vector of phenotypic observations for the traits that were assessed from each site (Tables 2-4); β is a vector of fixed effects, and this included the mean and missing values (mv) except where covariates were fitted (see results); u is the vector of random effects, which included replicates, blocks within replicates, tree (additive genetic effect-estimated using the relationship matrix derived from the pedigree file for trial trees and their ancestors), and specific combining ability (fitted by including the full-sib family identity in the model) terms; and e is a vector of random residuals. X and Z correspond to design matrices relating the observations in y to the fixed and random effects in β and u, respectively. The joint distribution of the random terms was assumed to be multivariate normal, with means and (co)variances defined as: where 0 is a null matrix and G and R are (co)variance matrices for effects in U and e, respectively [52][53][54]. Where rough bark height was fitted in the model, all the trees without rough bark were fitted as missing values for the trait being modelled.

Spatial Analyses
The presence of spatial heterogeneity in data may cause inaccurate estimation of genetic parameters [55]. There are diverse approaches for accounting for spatial effects [56]. However, to detect spatial effects in this study, a spatial term was introduced in the linear mixed models by fitting a two-dimensional spatial term [52,56,57]. To set the spatial term, every tree was uniquely identified by a row and column position within each trial, setting the absent, dead, and filler trees to missing values [52]. The missing values (mv) were included as a fixed factor in the models [53]. The spatial term was then partitioned into spatially correlated (ξ) and uncorrelated (η) residuals. The spatially correlated error (ξ) was modelled using a first-order separable autoregressive model in the row and column directions [44,51,58]. However, in addition to the two-dimensional separable first-order autoregressive spatial model, an independent residual (nugget-ψI150) was also added as a random term with the following form: where σ 2 is the spatial variance, ⊗ is the Kronecker product, and (p) is a first-order autoregressive correlation matrix with autocorrelation p for columns (col) and rows (row) [51].
To test the improvement of the models after fitting the spatial autocorrelation structure to the residuals, the spatial model was compared with the corresponding reduced model that did not fit a

Linear Models
All the univariate and bivariate linear mixed models were run in ASReml-R [51]. First, univariate models were fitted to generate variance components that were used to test the presence of additive genetic variation. The general linear mixed model is represented as: where y is the vector of phenotypic observations for the traits that were assessed from each site (Tables 2-4); β is a vector of fixed effects, and this included the mean and missing values (mv) except where covariates were fitted (see results); u is the vector of random effects, which included replicates, blocks within replicates, tree (additive genetic effect-estimated using the relationship matrix derived from the pedigree file for trial trees and their ancestors), and specific combining ability (fitted by including the full-sib family identity in the model) terms; and e is a vector of random residuals. X and Z correspond to design matrices relating the observations in y to the fixed and random effects in β and u, respectively. The joint distribution of the random terms was assumed to be multivariate normal, with means and (co)variances defined as: where 0 is a null matrix and G and R are (co)variance matrices for effects in u and e, respectively [52][53][54].
Where rough bark height was fitted in the model, all the trees without rough bark were fitted as missing values for the trait being modelled.

Spatial Analyses
The presence of spatial heterogeneity in data may cause inaccurate estimation of genetic parameters [55]. There are diverse approaches for accounting for spatial effects [56]. However, to detect spatial effects in this study, a spatial term was introduced in the linear mixed models by fitting a two-dimensional spatial term [52,56,57]. To set the spatial term, every tree was uniquely identified by a row and column position within each trial, setting the absent, dead, and filler trees to missing values [52]. The missing values (mv) were included as a fixed factor in the models [53]. The spatial term was then partitioned into spatially correlated (ξ) and uncorrelated (η) residuals. The spatially correlated error (ξ) was modelled using a first-order separable autoregressive model in the row and column directions [44,51,58]. However, in addition to the two-dimensional separable first-order autoregressive spatial model, an independent residual (nugget-ψI 150 ) was also added as a random term with the following form: where σ 2 is the spatial variance, ⊗ is the Kronecker product, and (p) is a first-order autoregressive correlation matrix with autocorrelation p for columns (col) and rows (row) [51].
To test the improvement of the models after fitting the spatial autocorrelation structure to the residuals, the spatial model was compared with the corresponding reduced model that did not fit a spatially correlated residual using a two-tailed likelihood ratio test (LRT) with 3 degrees of freedom [57].
For a visual inspection of the presence of small-scale spatial effects, semi-variograms of the residuals from the above spatial models were generated in R v. 3.6.0 [59] using the sp and gstat packages. The variograms were generated for selected scores of bark stripping and height, but the importance of spatial models was tested for all the variables that were assessed (see results). Flat variograms are expected for randomly distributed data. If spatial dependence is present, semivariance will be small at short distances, will increase at intermediate distances, and will reach an asymptote at longer distances [55]. The direction of strongest spatial correlation in each site was used to generate the variogram although there was spatial autocorrelation in all directions (results not shown). Of the five common spatial models (linear, exponential, spherical, Gaussian, and Matern), the Matern model was selected as the most suitable, suggesting the presence of irregular spatial patterns in both height and bark stripping [60].

Estimation of Additive Genetic Variation and Heritability within Sites
Because of spatial effects (see results), spatial models were used to obtain variance components to enable testing of the presence of additive genetic variation and for specific combining ability (SCA) for each trait that was assessed [51]. To test whether the additive genetic or SCA variation were greater than zero, the reduced models were compared with respective full models using one-tailed log likelihood ratio tests (LRT) with one degree of freedom [51]. For open-pollinated Payanna and Beulah trials where the SCA term was not fitted in the model, the reduced models excluded the additive genetic variation. For the control-pollinated Wilmot trial, where the additive and SCA terms were in the model, the random terms were sequentially dropped from the model, starting with the SCA term and then the additive genetic term.
For the binary trait, rough bark, a generalized linear mixed model with a binomial distribution and logistic link function was used to estimate variance components [51]. The importance of the random factors in the binomial models were evaluated using the Akaike information criterion [61].
For all models, individual narrow-sense heritability (h 2 ) was estimated from univariate spatial models, that included all the design terms as described above as the additive genetic variance (σ a 2 ) divided by the sum of σ a 2 and spatially independent error variance (σ e 2 ) [51] as below: The SCA term was excluded in the final models of the Wilmot trial because it was not significant for all traits (see results; Table 4). For the binary trait, the residual error variance on the underlying logistic scale, which is 3.28987 [51], was used in the heritability estimations. Estimates of the associated standard error for the estimated heritability were obtained from the average information matrix, using a Taylor series approximation [51]. The presence of additive genetic variation for survival was not tested due to the high survival in all the three trials.
Bivariate models were then fitted to determine the correlation of traits between trials (Type B) and within trials (Type A). Equation (1) was modified into a general bivariate model below: where y 1 and y 2 = vectors of measurements; X 1 and X 2 = incidence matrices of 0's and 1's, relating measurements and fixed effects; β 1 and β 2 = unknown fixed quantities including the underlying means for each trait; Z 1 and Z 2 = incidence matrices relating measurements and random effects; u 1 and u 2 = vectors of the unknown additive genetic effects; and e 1 and e 2 = vectors of unknown random errors.

Type B Genetic Correlations
In evaluating the genotype by environment (G×E) interaction, type B genetic correlations are used to assess the relative performance of traits measured in different environments [62]. The across-site additive genetic correlation of bark stripping, height, DBH and bark traits was assessed for Payanna and Beulah. This was carried out by treating the same trait at the two sites as different traits in a bivariate model [63]. The fixed term was only the mean while the random terms comprised the additive genetic (i.e., tree) term and site-specific replicates and incomplete blocks terms. The unstructured design matrix fitted allowed covariation between additive genetic variation, but a diagonal matrix was fitted for replicates and incomplete block levels. Spatial models were not used at this stage. The additive genetic correlation (r a ) was estimated as: where cov a (x,y) is the additive genetic covariance between traits x and y, σ 2 ax is the additive genetic variance components for trait x, and σ 2 ay is the additive genetic variance components for trait y. Standard errors of the genetic correlations were estimated as mentioned above [51].
To test whether the genetic correlations were significantly different from zero, a full model was compared with the respective reduced model that had the additive covariances fixed to zero using a two-tailed LRT with one degree of freedom. To test whether genetic correlations were equal to one (perfect correlation), a full model was compared with the respective model that had the additive covariances fixed to one using a one-tailed likelihood ratio test (LRT) with one degree of freedom.

Phenotypic and Type A Genetic Correlations
Within sites (Type A), additive genetic (r a ) and Pearson's phenotypic correlations among traits were estimated directly from non-spatial bivariate models as defined above (Equation (4)) with the mean as the fixed term and the design (replicates and incomplete blocks) as well as the additive genetic effects as random terms. The SCA term was not included in the model used for the Wilmot estimates. The models allowed covariation at the additive genetic variation, residual, replicates, and incomplete block levels. For binomial traits, the genetic correlations were estimated on the observed binary scale [64]. To test whether the genetic correlations were significantly different from zero, a full model was compared with the model where additive covariance was fixed to zero using a two-tailed LRT with one degree of freedom [51]. Testing that the Pearson's phenotypic correlations were not equal to zero was performed in R using the function cor.test.
To test if the measured physical traits fully explained the additive genetic variation in bark stripping, linear models (Equation (1)) were re-run with the traits that significantly correlated with bark stripping included as covariates (see results). As previously mentioned, the significance of the additive genetic variation in bark stripping from zero was tested with one-tailed LRT. The random term included the replicates, blocks within replicates and the additive genetic variation.
To test the performance consequences of bark stripping at Wilmot, genetic correlations of bark stripping assessed at year 2 and year 3 with the diameter at year 5 and change in height between consecutive years were estimated. This height change was assessed as the difference between height measurements taken in the different years, where ∆Ht2−3 = height year 3 minus height year 2 and ∆Ht3−5 = height year 5 minus height year 3.

Estimation of Genetic Gain
To determine the family rankings for the estimation of genetic gain for a single backward selection event, family best linear unbiased predictors (BLUPs) were estimated for all trials using univariate family models with the spatial term. The additive tree term was not included in the model at this stage, and thus the differences among families in these new models refers to the among family variance directly, which in the case of Wilmot includes the additive genetic and SCA effects estimated in the initial model. Genetic gain was estimated at the family level since the parents of the families are already available in seed orchards and can be used to generate new seed for selected families for deployment.
The predicted percentage genetic gain (reduction in bark stripping) from the backward selection of families was calculated as the average BLUPs of the selected less susceptible families for each trial divided by the population mean, multiplied by 100 [53]. Genetic gain was predicted for selection of a different proportion of families at each of the sites. Predictions were based on assessments at year 5 (2016) for Payanna and Beulah and the year 2 assessment (2017) for Wilmot that exhibited the highest heritability estimates.

Differences between Sites in Bark Stripping and Associated Traits
The percentage of trees exhibiting bark stripping was variable across sites, with 95%, 70%, and 52% of the trees experiencing some level of back stripping at Beulah (n = 2002), Payanna (n = 2668) and Wilmot (n = 1372), respectively ( Figure 3). Consistently, the amount of bark removed at Beulah (Table 2) was higher than the other two sites (Tables 3 and 4). Within the two older trials, the mean bark stripping was greater in the second assessment (Tables 2 and 3). However, in the younger Wilmot trial lower levels were observed in the second assessment, possibly due to fast wound recovery in young trees. Beulah, which had the highest levels of bark stripping, had trees that were shorter at age 4 years (x = 391 ± 89 cm) than the sister trial at Payanna (x = 573 ± 94 cm) (Tables 2 and 3). The DBH of the trees in Beulah was also lower than that of Payanna. The predicted percentage genetic gain (reduction in bark stripping) from the backward selection of families was calculated as the average BLUPs of the selected less susceptible families for each trial divided by the population mean, multiplied by 100 [53]. Genetic gain was predicted for selection of a different proportion of families at each of the sites. Predictions were based on assessments at year 5 (2016) for Payanna and Beulah and the year 2 assessment (2017) for Wilmot that exhibited the highest heritability estimates.

Differences between Sites in Bark Stripping and Associated Traits
The percentage of trees exhibiting bark stripping was variable across sites, with 95%, 70%, and 52% of the trees experiencing some level of back stripping at Beulah (n = 2002), Payanna (n = 2668) and Wilmot (n = 1372), respectively ( Figure 3). Consistently, the amount of bark removed at Beulah ( Table 2) was higher than the other two sites (Tables 3 and 4). Within the two older trials, the mean bark stripping was greater in the second assessment (Tables 2 and 3). However, in the younger Wilmot trial lower levels were observed in the second assessment, possibly due to fast wound recovery in young trees. Beulah, which had the highest levels of bark stripping, had trees that were shorter at age 4 years ( ̅ = 391 ± 89 cm) than the sister trial at Payanna ( ̅ = 573 ± 94 cm) (Tables 2 and  3). The DBH of the trees in Beulah was also lower than that of Payanna.
The bark traits assessed were different between the older (Beulah and Payanna) and younger (Wilmot) trials. The young Wilmot trees had no rough bark and the bark was too thin to be assessed even at the time of the 2nd assessment at age of 3 years. At 4 years of age, 1% of the trees in Payanna had rough bark compared with 5.3% at Beulah, but by 5 years of age 76% and 42.6%, respectively, of the trees had rough bark. Of the trees with rough bark at Payanna and Beulah, the mean height of rough bark on Payanna trees ( ̅ = 93.5 ± 60.4 cm) was higher than at Beulah ( ̅ = 47.9 ± 64.3 cm) (Tables 3 and 4). These differences in bark development at age 5 years may reflect different growth rates, with the faster growing trial at Payanna having a greater rough bark development.  The bark traits assessed were different between the older (Beulah and Payanna) and younger (Wilmot) trials. The young Wilmot trees had no rough bark and the bark was too thin to be assessed even at the time of the 2nd assessment at age of 3 years. At 4 years of age, 1% of the trees in Payanna had rough bark compared with 5.3% at Beulah, but by 5 years of age 76% and 42.6%, respectively, of the trees had rough bark. Of the trees with rough bark at Payanna and Beulah, the mean height of rough bark on Payanna trees (x = 93.5 ± 60.4 cm) was higher than at Beulah (x = 47.9 ± 64.3 cm) (Tables 3 and 4). These differences in bark development at age 5 years may reflect different growth rates, with the faster growing trial at Payanna having a greater rough bark development.

Spatial Effects
The variation in traits assessed was not randomly distributed within the blocks of the three field trials, as evidenced by the high significance of the spatial term of the various models (Tables 2-4). For bark stripping, more damage occurred at the edge of the blocks possibly due to cover for the marsupials provided by the windrows (Supplementary Figure S1). For purposes of illustration, Supplementary Figure S2 shows that plants in close proximity tended to be bark stripped at similar levels to those of plants that are spatially separated.

Additive Genetic Variation for Bark Stripping
There was significant (p < 0.001) additive genetic variation for bark stripping at all trials for at least one of the yearly measurements, and heritability estimates ranged between 0.06-0.14 (Tables 2-4). The lowest additive genetic variation and lowest heritability estimates were associated with the first bark stripping assessment at Payanna (p < 0.05; Table 3) and the second assessment in Beulah and the younger Wilmot trial (p < 0.05; Table 4). Within sites, the time of assessment was important in detecting significant additive genetic variation for bark stripping, and in part this appeared to reflect the intensity of damage. At Payanna, for example, the amount of bark stripping in the first assessment (age 4 years) was lower than in the following year (means of 8% versus 16%, respectively) (Tables 2 and 3), and exhibited less additive genetic variation and lower heritability, suggesting few trees being visited by the animals. However, bark stripping across years was highly additively genetically correlated at Beulah (r a = 0.78 ± 0.16), Payanna (r a = 0.91 ± 0.23) and Wilmot (r a = 0.99 ± 0.27), indicating either a confounding of the previous years' damage or that the family choice by the animals was consistent across years.
Most other traits showed significant additive genetic variation in one or all trials (Tables 2-4). There was significant additive genetic variation for height in all trials (p < 0.001). The trait with the highest heritability estimate in both of the older trials was rough bark height; for the cohort of trees that had rough bark and in the younger trial, bark strip height had the highest heritability. The SCA effects estimated from the full-sib family trial at Wilmot were insignificant for all traits, including bark stripping (Table 4).

Genetic × Environment Interaction
The additive genetic correlations for the same traits (Type B) assessed in the two sister trials-Beulah and Payanna-are shown in Table 5. There was no significant across-site genetic correlation for bark stripping at year 4 (r a = 0.23 ± 0.39), likely due to a low genetic signal with the low damage level observed at Payanna. However, with the cumulative and new damage at year 5, the heritability increased at Payanna and a significant and high positive across-site additive genetic correlation was observed for bark stripping (r a = 0.76 ± 0.25). This indicates a high correspondence of family ranks between the two trials at this age. Similarly, high correlations were obtained for the year 5 rough bark (r a = 0.74 ± 0.42) and the height at year 4 (r a = 0.91 ± 0.32). Table 2. Statistics for each trait assessed at Beulah. SD = standard deviation; h 2 (se) = narrow-sense heritability and its associated standard error estimated with univariate spatial models. The significance of the additive genetic variation and the spatial effect on the model are also presented. For the binary trait rough bark (marked *), + signifies that the parameter was important based on the Akaike Information Criteria. Genetic variation was not estimated (NE) for survival, since it was high for most families.

Trait
Year  Table 3. Statistics for each trait assessed at Payanna. SD = standard deviation; h 2 (se) = narrow-sense heritability and its associated standard error estimated with univariate spatial models. The significance of the additive genetic variation and the spatial effect on the model are also presented. For the binary trait rough bark (marked *), + signifies that the parameter was important based on the Akaike Information Criteria. Genetic variation was not estimated (NE) for survival, since it was high for most families.  Table 4. Statistics for each trait assessed at Wilmot. SD = standard deviation; h 2 (se) = narrow-sense heritability and its associated standard error estimated with univariate spatial models. The significance of the additive genetic variation and specific combining ability (SCA) and the spatial effect on the model are also presented. Genetic variation was not estimated (NE) for survival, since it was high for most families.  Table 5. Across-site (Type-B) additive genetic correlations (r a ) for each trait assessed in both Beulah and Payanna and the corresponding standard error (se). The two sister trials share 98 families and were the same age. The chi-square value (χ 2 ) and the p-value that the genetic correlation is equal to zero (P[r a = 0]) are indicated. When significantly different from zero (p < 0.05), the significance of the genetic correlation from 1 (P[r a = 1]) is also presented. A " + " indicates that the difference of the genetic correlation for the binary rough bark from 0 or 1 was important based on the AIC.

Traits Associated with Bark Stripping
The phenotypic (0.34 to 0.71) and genetic (0.78 to 1.00; Type A) correlations between the various measures or ages of bark stripping were generally highly significantly positively correlated within a trial (Tables 6-8). The phenotypic and genetic correlations of bark stripping with growth traits-height and stem diameter-varied from positive to negative (Tables 6-8). In general, bark stripping was significantly negatively correlated with stem diameter at the phenotypic level, and while this trend was still evident at the genetic level these correlations were not significant. The association with height was more variable. At the phenotypic level, the correlation was generally negative, but the genetic correlation tended to be negative in the older trials (Tables 6 and 7). The exception was the year 5 bark stripping at Beulah, which was significantly positively correlated with height (r a = 0.33 ± 0.34), possibly due to the tendency for rough bark height to negatively genetically correlate with height at this age (see below). In the younger Wilmot trial, the phenotypic correlation between bark stripping and height was evident at the onset of bark stripping (age 2 years), but this subsequently became significantly negative consistent, with an adverse effect of bark stripping. In contrast, the genetic correlations between bark stripping and height were generally positive, indicating that families that were initially faster growing had more bark stripping, although this was only significantly greater than zero in the case of 3-year bark stripping and 2-year height.
In the older trials, bark stripping was phenotypically significantly negatively correlated with year 5 bark thickness, the presence of rough bark, and the rough bark height, with one exception (Payanna bark stripping at age 4 years versus 5 year rough bark height, possibly because there were very few trees with rough bark at age 4) ( Tables 6 and 7). This negative association was also evident at the genetic level, with significant negative correlations of year 5 bark stripping with bark thickness (r a = −0.48 ± 0.22), rough bark (r a = −0.47 ± 0.19)" and rough bark height (r a = −0.37 ± 0.17) but only in Payanna ( Table 7), suggesting that these traits influence the level of bark stripping at the phenotypic and genetic levels. Compared to Beulah, we noted a rapid increase in the amount of rough bark between 4 and 5 years (77% of the trees in Payanna had rough bark in year 5 compared to 42% at Beulah), emphasizing importance of this trait in minimizing bark stripping, but this depends on its extent in the population. In the younger Wilmot trial, rough bark had not developed at the time of bark stripping. However, in this case, both the 2-and 3-year bark stripping and 2-year bark strip height were highly significantly positively correlated with stem access at the phenotypic (0.31 to 0.42) and genetic (0.87 to 1.00) levels (Table 8). Table 6. Phenotypic (above diagonal) and additive genetic Type A (below diagonal) correlations of traits assessed at the Beulah trial with the standard errors given in parentheses. The significance of the difference of the correlations from zero are indicated: *** = p < 0.001, ** = p < 0.01, * = p < 0.05. NA = not applicable, as only assessed on a specific cohort of individuals. The genetic correlation between bark stripping for year 2 and 3 and height or DBH for year 5 in Wilmot was negative, suggesting a tendency for bark stripping to reduce performance. A negative genetic correlation was detected between year 5 diameter and bark stripping for year 2 (r a = −0.84 ± 0.16) and year 3 (r a = −0.77 ± 0.22). Similarly, a nonsignificant negative correlation between the first height increment (∆Ht2-3) and bark stripping for year 2 (r a = −0.35 ± 0.27) and year 3 (r a = −0.26 ± 0.29) was detected, indicating a reduction in performance in the initial years of bark stripping, with recovery in the later years.
To test if the observed additive genetic variation in bark stripping could be fully explained by tree physical traits, linear models were run with traits that significantly genetically correlated with bark stripping (i.e., height, stem access, rough bark, and bark thickness) as covariates (Table 9). However, the additive genetic variation in bark stripping was still significant after accounting for this covariation at the phenotypic level. In the younger Wilmot trial, for example, after simultaneously including stem access as a covariate in the model for year 2 bark stripping, the additive genetic variation was still significant (LRT χ 2 = 15.6, p < 0.001) but heritability reduced from h 2 = 0.09 ± 0.03 to h 2 = 0.06 ± 0.02. Similarly, in Payanna, the additive genetic variation for bark stripping was significant after accounting for covariation with selected traits (Table 9). Table 9. Genetic estimates for bark stripping derived from models that accounted for covariation in tree physical traits in Payanna and Wilmot trials. No significant correlations between bark stripping and physical traits were detected in Beulah. Shown are the chi-square value (χ 2 ) and probability for the one-tailed likelihood ratio test (LRT) that the additive variance (V a ) in bark stripping is greater than zero and its narrow-sense heritability (h 2 ) and its associated standard error (se) estimated with univariate spatial models. The heritability estimates calculated without fitting covariates are given in Tables 2-4.

Estimation of Genetic Gain
Depending on the number of families selected, an up to 37% reduction in bark stripping was predicted to be achievable based on family selection using their best linear unbiased predictions (BLUPs) (Figure 4). However, selecting 20% of the least susceptible families in each trial resulted in an expected gain in the reduction in bark stripping of 7.9% at Beulah, 3.8% at Payanna, and 22.1% at Wilmot. The highest predicted genetic gain was from family selection in the younger, full-sib progeny at the Wilmot trial. The half sib-OP families at Payanna and Beulah gave relatively low genetic gain predictions.
Forests 2020, 11, x; doi: FOR PEER REVIEW www.mdpi.com/journal/forests predicted to be achievable based on family selection using their best linear unbiased predictions (BLUPs) (Figure 4). However, selecting 20% of the least susceptible families in each trial resulted in an expected gain in the reduction in bark stripping of 7.9% at Beulah, 3.8% at Payanna, and 22.1% at Wilmot. The highest predicted genetic gain was from family selection in the younger, full-sib progeny at the Wilmot trial. The half sib-OP families at Payanna and Beulah gave relatively low genetic gain predictions.

Discussion
Using three field-based genetic trials established with a large number of Pinus radiata families, we show the existence of additive genetic variation in susceptibility to bark stripping damage by marsupial herbivores, which, when well-expressed, appears relatively stable across sites. Plantation age, bark features, and the presence of needles or branches covering the stem were important correlates with the amount of bark removed from the trees by the marsupials. Importantly, these physical traits are under significant additive genetic control and can thus also be enhanced through selection. An up to 22.1% reduction in bark stripping can be achieved by selecting 20% of the less susceptible families, although the stability of the less susceptible families, when planted separately, needs further testing. This is the first study, to our knowledge, to demonstrate genetic gain through the selection of families less susceptible to bark stripping.
The estimated narrow-sense heritability of bark stripping (0.06 to 0.14) was within the lower range of what has been reported for the damage of conifer bark from insect herbivores (0.02-0.40) [12,13,23,[65][66][67]. Similarly, it is in the lower range of the narrow-sense heritabilities reported for resistance to needle pathogens (h 2 = 0.05 to 0.69) in Australian P. radiata plantations [18]. Estimates were also markedly lower than what has been reported for mammalian browsing on the needles of Pseudotsuga menziesii Mirb. (h 2 = 0.73) [29]. Since P. radiata has not coevolved with the marsupials and animals have generally had a relatively minor impact on the natural P. radiata stands in California, its place of origin, the presence of genetic variation could be a by-product of variation in resistance traits to other native biotic threats.
Within a particular site, the results indicated the stability of genetic variation in damage across different years depicted by the high positive genetic correlation of the damage between years. This may show that the families damaged in one year are consistently damaged in the subsequent year or this may be an effect of confounding damage from the previous year, as the subsequent scores included old damage. However, the plant traits that correlate with bark stripping are likely to change with age, as depicted by the marked increase in trees developing rough bark between 4 and 5 years in the older trials. In this case, the changes in rough bark are relatively consistent among the families for each year, but the stability of the genetic signal for bark stripping between periods of contrasting developmental stages (for example, at 2 years before rough bark development vs. 5 years) is unknown. In the different years, the genetic estimates are likely to be in part influenced by the amount of bark stripping. For most biotic stresses, the differentiation between susceptibility categories becomes clearer with the increase in damage intensity, since high damage levels reduce the possibility of random escape from damage [12,40,68,69]. However, cumulative damage estimates are contingent upon the survival of the trees and the rate of bark recovery after the initial damage. At the Wilmot site, for example, the cumulative bark stripping was lower possibly due to the faster recovery of the younger trees with low damage. The threshold damage to enable the differentiation of the genotypes is not known, since susceptible individuals may appear resistant if they are located in areas of low bark stripping intensity, and conversely the resistance in less susceptible genotypes can be overcome in regions of high-intensity damage. Therefore, to correctly identify less susceptible genotypes, selected families should be tested in areas of high bark stripping intensity to determine whether extreme browsing intensity does not lead to the browsing of resistant families.
When genetic variation is well expressed (age 5 years), the bark stripping in Beulah and Payanna was strongly genetically correlated, signaling low G×E and reflecting the fact that the genes influencing bark stripping and their expression are similar in different environments. Indeed, this is consistent with the negative genetic correlation between fifth year bark stripping and rough bark height in both trials, suggesting an underlying common mechanism. In this case, the relative resistance of different families was not dependent on the environment, so the families that received relatively less bark stripping in Payanna were less bark-stripped in Beulah. Studies have shown strong positive genetic correlations between trials for damage by insect herbivores in P. radiata [41] and in spruce [12,26,70], even when the pest was at a low density [26]. However, in the present case the high genetic correlations for bark stripping were only true for trials of the same age that had similar physical features to deter bark stripping (e.g., well-developed rough bark), and the extent to which G×E is affected when the families are at different development stages needs to be tested.
This study provides evidence that physical traits play a role in moderating bark stripping preferences in the field. The role of the physical traits appears to vary depending on the age of the trees. In older trees (>3 years), the physical traits of the bark (i.e., thick and rough bark) may offer protection against bark stripping. In the younger trees < 3 years, because the bark features are not well developed the presence of obstructive needles and low branches are a significant explanatory factor. In P. radiata, bark thickness is positively correlated with age, height, and DBH [71]. The negative association of bark stripping with rough bark development in P. radiata has been previously noted [6]. Studies in other conifers have also indicated the importance of rough and thick bark [72][73][74][75] and other physical features, such as obstructive branches on the stem [74], in deterring mammalian herbivores. Thick bark may be difficult to detach but has also been associated with a high density of resin canals and gritty-textured sclereids [76,77]. These traits often occur simultaneously, but the relative importance of each trait is not known. In contrast to the negative association observed for mammal damage, damage by insects has been positively correlated with rough bark in conifer species [78], which highlights that resistance strategies may be herbivore-specific. Rough bark supports wood-boring insects, as the fissures are a safe place for oviposition or protection from natural enemies and environmental extremes. It is possible that these organisms may distract bark stripping mammals which could further contribute to reducing susceptibility. The bark features were genetically variable, with possibly no genetic constraints to their selection indicated by the tendency for positive genetic correlations among these traits. Similarly, the bark features showed a low G×E interaction, and this may in part explain the low G×E for bark stripping.
Height was also genetically correlated with bark stripping, although the trend contrasted in the trials of different ages. In the young population (Wilmot), the herbivores appeared to be attracted to faster growing trees, in contrast to the two older plantations where faster growing trees were less damaged. A comparison of the relative size of trees damaged by insect and mammalian herbivores in other P. radiata studies [79] and in other conifer [13,67,70] and non-conifer [80] species suggests that trees most often damaged are those which are larger, dominant and growing most rapidly. Damage to larger trees has been correlated with high sap volume, phloem thickness, and sap sweetness [81]. Additionally, faster growing trees have been reported to invest less in defense [82]. In young trees, bigger trees possibly give more bark per unit stripping than smaller trees and could possibly produce more sugars. In contrast, the results showed that in older conifer plantations animals prefer smaller trees, possibly with less developed bark features, as evidenced by the positive genetic and phenotypic correlation between rough bark or bark thickness with height and DBH. It is also possible, however, that the smaller trees have remained thin as a result of consecutive bark stripping, depicted from the strong negative correlation of initial bark stripping with later DBH. The selectivity of smaller trees by bark stripping mammals is common in mature conifer plantations [72,74,83].
It has also been suggested that fast growing trees should be able to recover more quickly than slow growing trees [84], an aspect of tolerance. In this study, tolerance may be indicated by the positive correlation between early bark stripping and later height. Heritable variation in tolerance against mammalian and insect herbivory has been observed in other Pinus spp [85,86], where trees appear not to employ resistance strategies but counter the effects of the browsing through compensatory growth. For commercial P. radiata plantations intended for timber production, tolerance will be less desirable, since bark stripping exposes tissues to fungal attack with subsequent rotting, which may reduce timber quality [30]. In contrast with height, there was a negative genetic correlation between early bark stripping and later diameter, suggesting that heavily damaged trees were putting on less volume than the less damaged families. However, the presence of residual variation, which could not be explained by physical or growth traits, suggests the possible involvement of chemical features in driving differences in bark stripping. This is yet to be tested.
While the capturing of genetic gain in reduced bark stripping could be achieved by cloning selected trees based on their breeding or deployment values, this would require rejuvenation, vegetative propagation, and subsequent clonal testing [87], which will have a long time lag. Given that the families are derived from seed orchard parents, the rapid deployment of families of low susceptibility to bark stripping is an obvious path to rapidly capture genetic gain from the current trials. However, the predicted genetic gains from such backward family selection assumes that the family ranking for bark stripping do not change in different environments, as this study suggests. In this case, the genetic gain for reducing susceptibility to bark stripping is estimated to be between 3.8% and 22.1% when 20% of the most resistant families are selected for deployment. This is in the genetic gain range that has been estimated for herbivory in conifers [29]. In Pseudotsuga menziesii, based on selecting 10% of the most resistant families against deer browsing, genetic gain was estimated at 11% [29]. Similarly, against the pine weevil in Picea abies H. Karst. genetic gain varied between 8% and 50% [12]. However, the genetic gain will be influenced by the proportion of the less susceptible families selected and applies to current field setting. Whether or not the less susceptible families retain their rank when grown separately needs to be tested. According to the associational resistance theory, resistant plants may appear more resistant in the presence of a more susceptible plant that can easily be eaten by the animals [88]. If the marsupials use the bark as a supplementary rather than an alternative feed, as suggested by Smith, Ratkowsky [49], the damage in monocultures may persist. Alternatively, deploying more resistant P. radiata plants may shift the browsing pressure onto other more susceptible non-plantation species in the coupe [49], which could have additional positive effects by reducing interspecific plant competition with the P. radiata [89].

Conclusions
This paper illustrates that variation in bark stripping is under low but significant additive genetic control and, when well expressed, the genetic signal appears to be relatively stable across different environments. This provides an opportunity for selection for reduced susceptibility with potential genetic gains for deployment and breeding. Based on phenotypic and genetic correlations, several physical plant traits were identified as being likely to contribute to the variation in bark stripping. Initial variation in damage may be affected directly by plant size and accessibility, but later factors such as bark thickness and bark texture become important. The presence of unexplained genetic variation after accounting for these physical factors suggests that other characteristics, such as endogenous chemicals, may also be involved in explaining the variation in the bark stripping of P. radiata.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4907/11/12/1356/s1: Figure S1: Heat maps of bark stripping scores from the three field sites based on assessments for year 5 (2016) at Beulah and Payanna and for year 2 (2017) at Wilmot, Figure S2: Semivariograms of the residuals generated from the spatial models showing the small-scale spatial variation in bark stripping (above) and height (below) for (a) Beulah, (b) Payanna. and (c) Wilmot after adjusting for random genetic, replicate. and incomplete blocks within replicate effects using mixed models. If spatial dependence is present, semivariance will be small at short distances, will increase at intermediate distances, and will reach an asymptote at longer distances. Estimates are based on assessments at year 5 (2016) at Beulah and Payanna and at year 2 (2017) for Wilmot, except for the height variogram for Payanna that used the year 4 height, since only one height measurement was taken for this trial.