Introduction
Genetically-based resistance to insect pests is well documented for many agricultural crop species (Maxwell et al. Reference Maxwell, Jenkins and Parrott1972). However, despite the large number of forest tree breeding programmes, little is known about the heritability of resistance against insects in coniferous trees. Forest geneticists have evaluated the heritability of many traits of conifers but most of this research was focussed on traits associated with growth and wood properties. Heritability of insect resistance in conifers has only been reported for a few species of herbivorous insects, including moths (Lepidoptera) (Quiring et al. Reference Quiring, Turgeon, Simpson and Smith1991), sawflies (Hymenoptera) (Leinekugel le Cocq et al. Reference Leinekugel le Cocq, Quiring, Verrez and Park2005), aphids (Hemiptera: Aphididae) (Jensen et al. Reference Jensen, Harding and Roulund1997), adelgids (Hemiptera: Adelgidae) (Nielsen et al. Reference Nielsen, Kirkeby-Thomsen and Roulund2002; Bains et al. Reference Bains, Isik, Strong, Jaquish, McLean and El-Kassaby2009), and weevils (Coleoptera: Curculionidae) (King et al. Reference King, Yanchuk, Kiss and Alfaro1997; Zas et al. Reference Zas, Sampedro, Prada and Fernández-López2005; Verrez et al. Reference Verrez, Quiring, Leinekugel Le Cocq, Adams and Sung Park2010).
Balsam fir, Abies balsamea (Linnaeus) Miller (Pinaceae), is a common tree species in the forests through northeastern North America. It is an ecologically and economically important species making up a large proportion of the boreal and Acadian forest in this area (Bakuzis and Hansen Reference Bakuzis and Hansen1965). This species is harvested mainly for lumber and pulpwood but is also used for Christmas trees and holiday greenery. Reforestation of harvested balsam fir stands is achieved primarily through silvicultural techniques used to promote natural regeneration, which in eastern North America yields between 6150 and 48 025 stems/ha (Pothier Reference Pothier2002). Because it regenerates prolifically, balsam fir has not been included in tree improvement programmes and few seedlings are cultivated for use in reforestation (Bakuzis and Hansen Reference Bakuzis and Hansen1965). However, Christmas trees are grown in intensively managed balsam fir stands with over 95% coming from plantations (Chastagner and Benson Reference Chastagner and Benson2000) and thus breeding strategies aimed at reducing pest damage would be beneficial to Christmas tree growers.
The balsam twig aphid, Mindarus abietinus Koch (Hemiptera: Aphididae), feeds on Abies Miller in North America and is a major pest in Christmas tree plantations of balsam fir and Fraser fir, Abies fraseri (Pursh) Poir, (Fondren and McCullough Reference Fondren and McCullough2003). Twig aphids feed on current-year shoots as they develop, causing needles to twist and distort, leading to severe aesthetic damage and reduced growth (Saunders Reference Saunders1969). Christmas tree growers frequently use pesticides to reduce aphid populations and mitigate damage (Fondren and McCullough Reference Fondren and McCullough2003). A significant proportion of variation in twig-aphid damage has been attributed to provenance for both balsam fir and white fir, Abies concolor (Gordon and Glendinning) Lindley ex Hildebrand (DeHayes Reference DeHayes1981; Ferrell Reference Ferrell1989), suggesting that there may be genetically-based resistance in firs to this pest.
Here we describe results of a four-year field study, carried out to determine: (i) the heritability of resistance of balsam fir to damage caused by the balsam fir twig aphid; and (ii) whether genetically-based resistance is consistently expressed by clones among years.
Materials and methods
Clonal seed orchards
We evaluated the susceptibility of balsam fir genotypes to the balsam fir twig aphid in two clonal seed orchards located in Atlantic Canada. Trees in each orchard were 5–10 m tall and reproductively mature. The orchard in Kingsclear, New Brunswick, Canada (45.9661°N, 66.7999°W) was established in 1992 by the New Brunswick Tree Improvement Council. It contained 24 clones with 6–10 ramets (x̄=8) per clone planted randomly throughout the orchard. Tree spacing was 3 m within rows and 6 m between rows. Clones were grafted from selected parent trees located throughout New Brunswick. Parent trees were plus-trees; a plus-tree is one that has been selected based on a subjective judgment of its phenotype but has not been tested to determine if desirable traits are heritable. Orchards were established for Christmas tree breeding purposes and initial parent tree selection was based on a number of traits, including growth rate, foliage arrangement, colour, needle retention, and branch form. Damage from insect pests was not directly considered during selection, although trees showing signs of heavy damage would not have been selected. Parent trees were selected from different stands and are not likely to be closely related.
The second orchard located in Debert, Nova Scotia (45.4172°N, 63.4117°W) was established between 1994 and 1997 by the Nova Scotia Department of Natural Resources. This orchard contained 187 clones, planted randomly with rows separated by 5 m and trees within rows by 2.5 m. Clones were from selected parent trees (plus-trees) within Nova Scotia; selection criteria were similar to that mentioned above for New Brunswick. Twenty-five clones were used in this site; clone selection was based on the number of ramets available and only clones that had ≥ 9 separate ramets were used.
Damage assessment
Aphid damage was assessed in the fall or early winter each year from 2009 to 2012. All trees were sampled at Kingsclear whereas 9–10 ramets of selected clones were sampled at Debert. Each year, one mid-crown branch from each ramet was cut from the southwest side of the tree. Branches were cut just below the fifth node along the main axis of the branch in 2009, and just below the third node in the subsequent three years. The current-year foliage of each branch was inspected and the percentage of shoots with distinctive needle curling, caused by aphid feeding, was recorded. Aphid damage at Debert was negligible in 2011 and 2012 and thus data for these two years were not included in analyses.
Statistical analysis
Traditionally most analyses of heritability use random (or mixed) effects linear models to estimate variance components within or between clones. While this is appropriate for normally distributed continuous traits, our data are proportions and are bound by zero and one. While many studies transform these data to perform standard tests, we chose statistical tests that would account for a binomial distribution. Our statistical analyses were carried out using generalised linear mixed models with a Bayesian estimation framework. We used the MCMCglmm package (Hadfield Reference Hadfield2010) in R version 3.0.3 (R Development Core Team 2014), which uses Markov chain Monte Carlo (MCMC) methods to estimate distributions of variance components. This method allows both a binomial response (i.e., the number of non-damaged shoots was considered a “success” and the number of damaged shoots was considered a “fail”) and the inclusion of both fixed (year) and random (clone and tree) effects. They also account for overdispersion (using an additive model) and allowed for specification of correlation matrices for the variance components (Hadfield Reference Hadfield2010). Overdispersion can occur if variance exceeds that predicted by the model. In additive models an additional random term is included into the model to capture any overdispersion (Browne et al. Reference Browne, Subramanian, Jones and Goldstein2005).
To determine the most appropriate model for our data we first modelled all the variables (fixed effect of year and random effects of clone, tree, and the clone×year interaction) in a single full model. We then compared sequentially reduced models (each with one term removed) to our full model using deviance information criterion values and examined variance distribution plots to assess overall model performance. From this we determined that a final model without the effect of tree or the clone×year interaction was the most appropriate. All of our models were fit with uninformative priors for both fixed and random effects and ran with a 10 000 burn-in period, followed by 100 000 iterations and a thinning rate of 100.
In each year we had different levels of aphid damage within our orchards. Previously, both Strong et al. (Reference Strong, Larsson and Gullberg1993) and Leinekugel le Cocq et al. (Reference Leinekugel le Cocq, Quiring, Verrez and Park2005) reported a positive relationship between heritability of host plant resistance and insect density. This suggests that at both very high and low densities of insects, genetically-based resistance may not be detectable, as variance among trees would be reduced. We wanted to estimate heritability for each year individually so we needed separate estimates of variance for each year. By default in MCMCglmm random variables are assumed to have a structured identity matrix (with ones on the diagonal and zeros on the offsets), this assumes independence among the levels but only one variance component is estimated (Hadfield Reference Hadfield2010). To estimate variance components for each year individually we changed the variance-covariance matrices for the random (clone) and residual variance in our model. Clone variance was expected to vary by year but also co-vary among years so an unstructured variance-covariance matrix was used, which results in one estimate of clone variance for each year and provides estimates of covariance between years. Residual variance was also estimated for each year individually using a covariance structure fixed to zero, as we had no expectation that residual variance would be correlated among years.
MCMCglmm models produce posterior distributions of estimated variance components (as opposed to point estimates). These variance distributions were used to estimate broad-sense heritability, the proportion of total genetic variation to total phenotypic variation (Vg/Vp), at both the individual
$({\rm H}_{i}^{2} )$
and clone-mean
$({\rm H}_{c}^{2} )$
level for each year for both the New Brunswick and Nova Scotia orchards following White et al. (Reference White, Adams and Neale2007):
where
$\sigma _{c}^{2} $
,
$\sigma _{\epsilon }^{2} $
, and r are, respectively, the estimated clone variance, estimated random error, and harmonic mean number of ramets per clone. Random error of a GLMM is not directly given but can be taken as:
where
$ {\rm \sigma }_{e}^{2} $
is the additive dispersion component and
${\rm \pi }^{2} /3$
is the distribution specific variation for a binomial model (Nakagawa and Schielzeth Reference Nakagawa and Schielzeth2010). Because we used variance distributions and not point estimates we created a posterior distribution of heritability estimates. The mode of the distribution was used to estimate heritability and 95% credible intervals were calculated.
To determine if clone performance (i.e., low or high aphid damage) is consistent over time we examined the correlation of damage levels over time. Pearson’s phenotypic correlations based on clone means, and genetic correlations of damage between all possible pair-wise combinations of years were estimated. Since aphid damage was measured on the same individual trees we used a type A (or “age-age”) genetic correlation. We treated aphid damage in each year as a separate trait, as in Falconer and Mackay (Reference Falconer and Mackay1996):
where x and y are two traits (twig-aphid damage in two separate years), cov
xy
the clonal covariance between the traits, and
${\rm \sigma }_{{cx}}^{2} $
and
$${\rm \sigma }_{{cy}}^{2} $$
are the clonal variance components for each trait. Estimates of variance and covariance were taken as the mode of the corresponding posterior distributions given by our statistical mode. Standard errors for the genetic correlations were not calculated since sampling errors are large and difficult to specify (Falconer and Mackay Reference Falconer and Mackay1996).
Results
In New Brunswick the mean percentage of shoots damaged by twig aphid decreased during our four-year study from 27.1% in 2009 to 6.3% in 2012 (Fig. 1). Overall, estimates of individual heritability were very low and varied between 0.03 (creditable interval (CI): 0.02–0.07) and 0.06 (CI: 0.03–0.16); clone explained little of the observed phenotypic variation within each year. Estimates of clone-mean heritability were moderate, and varied between 0.20 (CI: 0.08–0.48) and 0.50 (CI: 0.23–0.64) (Fig. 1). The 95% creditable intervals increased as damage within the orchard decreased, suggesting estimates of heritability are more robust as damage increases.
Relationship between (A) clone mean and (B) individual broad-sense heritability estimates (±95% credible intervals) and mean percentage of shoots damaged per tree each year by the balsam twig aphid in one clonal balsam fir orchard, located in New Brunswick (∙) and one in Nova Scotia (Δ).

In Nova Scotia the mean percentage of shoots damaged was low in 2009 and 2010 (8.0% and 1.3% respectively) and estimates of heritability were comparable to what was found in the New Brunswick orchard (Fig. 1). In 2010, posterior distribution plots showed that variance estimates for the effect of clone were stuck on zero and the chain did not converge. As we could not get reliable estimates of variance components for this year, we considered them to be effectively zero, and consequently heritability estimates were found to be near zero for Nova Scotia in 2010.
In New Brunswick, both pair-wise phenotypic and genetic correlations between the mean damage per clone by the twig aphid in one versus another year were positive, although not always significant (Fig. 2). In Nova Scotia, there was no significant relationship between mean clone damage in 2009 and 2010 (r p=−0.089, P=0.672).
Pair-wise correlations of mean clone damage from the balsam twig aphid between years in a New Brunswick balsam fir orchard. Both phenotypic (r p) and genetic (r g) correlations are reported where appropriate. Significant correlations (P<0.05) are indicated by an asterisk for r p; significance was not calculated for r g (see text for details).

Discussion
There was a moderate genetic basis in balsam fir for resistance to damage from balsam twig aphid. Individual estimates of heritability were low, suggesting that susceptibility to damage is not consistently expressed among ramets within a clone. Low individual heritability suggests that resistance to aphid damage is a phenotypically plastic trait. However, it is also possible that patchy distribution of aphids obscured reliable estimates of heritability especially when damage levels were low. Estimates of clone mean heritability were moderate. These estimates are naturally higher than those of individual heritability since within clone variance (error variance) is reduced while the variance between clones is unchanged. Moderate estimates of clone mean heritability suggest that despite the variation among individual ramets, the mean phenotypic value of resistance among clones is moderately influenced by genetics and that selection for resistance to aphid damage could result in genetic gains.
In the New Brunswick orchard, where damage by the twig aphid was highest, correlations of mean clone damage among years were positive and generally significant, which indicates that the resistance ranking of the different clones remained similar among years. Twenty-five per cent of clones consistently ranked in the top 50% of clones during each of the four years, further suggesting that selection for increased resistance is possible.
Low to moderate estimates of heritability may be a reflection of the complexity of herbivore-host interactions. Resistance to herbivores is a multivariate trait that is often influenced by many constitutive and induced defence mechanisms (see Hanover Reference Hanover1975; Larsson Reference Larsson2002). If multiple traits function together to influence tree resistance to aphid damage, each trait would presumably show some degree of phenotypic plasticity and each will be influenced by environmental (non-genetic) factors. Price and Schluter (Reference Price and Schluter1991) theoretically demonstrated that life-history traits, presumably influenced by many individual metric traits (e.g., morphological, behavioural, and physiological), are expected to show lower levels of heritability when compared to metric traits. They suggest that environmental variance accumulates from each metric trait involved, increasing environmental variance of life-history traits and thus reducing heritability estimates. Classical theory also suggests that fitness related traits may be under strong selective pressure with stabilising selection reducing the genetic variation among individuals and thus decreasing genetic variance and heritability estimates (see Mousseau and Roff Reference Mousseau and Roff1987).
We do not know which traits confer resistance of balsam fir to twig-aphid damage. Variation in tree susceptibility could be due to variations in the preference of the twig aphid for susceptible clones, similar to what has been reported for aphids on melon (Cucumis melo Linnaeus; Cucurbitaceae) (Kennedy and Kishaba Reference Kennedy and Kishaba1977), alfalfa (Medicago sativa Linnaeus; Fabaceae), and sweetclover (Melilotus Linnaeus; Fabaceae) (Kishaba and Manglitz Reference Kishaba and Manglitz1965), or to variation in aphid performance on different clones, as seen in cottonwood hybrids (Larson and Whitham Reference Larson and Whitham1997) and tomato plants (Kaloshian et al. Reference Kaloshian, Kinsey, Ullman and Williamson1997). Alternatively, variations in the susceptibility of balsam fir to the twig aphid may be due to genetically-based variations in tolerance of the host tree. Nettleton and Hain (Reference Nettleton and Hain1982) found that Fraser fir could outgrow up to 55% of twig-aphid damage. Similarly, large variations in the amount of damage found between groups of related (half-sibling) Nordmann fir (Abies nordmanniana (Steven) Spach (Pinaceae)) exposed to similar densities of silver fir woolly adelgid Dreyfusia nordmannianae (Eckstein) (Hemiptera: Adelgidae), which causes needle curl damage similar to that of twig aphid, suggests that there is genetically-based variation in tolerance of the host tree to this pest (Nielsen et al. Reference Nielsen, Kirkeby-Thomsen and Roulund2002).
If resistance of balsam fir to twig-aphid damage were in part due to tree tolerance, this would be beneficial to breeding programmes, as tolerance does not put selective pressure on pests to adapt and overcome host resistance (Rausher Reference Rausher2001). If balsam fir is able to outgrow twig-aphid damage our heritability values may be underestimated, because trees used in this study were previously selected as “plus-trees”. Parent trees were largely chosen for desirable growth characteristics (i.e., internodal and shoot length) and selecting trees for growth characteristics could have reduced the variation in resistance to aphid damage if resistance is correlated with growth related traits. Further research exploring the mechanisms of balsam fir resistance to aphid damage and the relative influences of insect preference, performance, and tree tolerance should be explored.
This study was conducted on older mature trees and heritability estimates can be influenced by organism age. Age dependent variation in genetic expression and heritability has been observed in animals (e.g., Charmantier et al. Reference Charmantier, Perrins, McCleery and Sheldon2006; Serbezov et al. Reference Serbezov, Bernatchez, Olsen and Vøllestad2010) and plants (e.g., Boege and Marquis Reference Boege and Marquis2005). Serbezov et al. (Reference Serbezov, Bernatchez, Olsen and Vøllestad2010) suggests that increases in variation with age may be attributed to differences in environmental factors that accumulate over time, resulting in greater phenotypic variation of traits. In trees these environmental effects could be due to a number of differing conditions including (but not limited to) microclimate, herbivore damage, pathogen/fungus attack, mechanical damage (e.g., ice or wind), and resource competition. Accumulated differences in environmental effects over time may reduce estimates of trait heritability in older trees. Also, resource allocation to plant defences can change over time during plant development (Boege and Marquis Reference Boege and Marquis2005). As trees used in this study were much older than marketable Christmas trees, further work will be needed to verify that the estimates of heritability of balsam fir to twig-aphid damage observed in this study are similar to that of younger trees.
The present study was replicated temporally but not spatially, as clones planted in New Brunswick were not the same as the ones planted in Nova Scotia. Environmental conditions can influence host plant phenotype and pest populations, and thus trees found to be resistant in one location may not be resistant in another if resistance is only expressed under certain conditions (Fritz Reference Fritz1990). Significant genotype×environment interactions in the resistance of conifers to other pests has been previously reported (Burdon Reference Burdon1977; Quiring et al. Reference Quiring, Turgeon, Simpson and Smith1991; Alfaro et al. Reference Alfaro, Borden, King, Tomlin, McIntosh and Bohlmann2002) and should be investigated in any tree-breeding programme.
In summary, our study provides evidence that incorporating genetically-based resistance into pest management strategies may be possible for reducing twig aphid-damage within balsam-fir Christmas tree stands. Further work is needed to determine the mechanisms of resistance to aphid damage and how these traits are expressed in younger (market-sized) trees. Also, future work should determine if selection for resistance to the twig aphid is compatible with selection for other commercially desirable traits, such as branch form, growth, needle retention, and resistance to other pests.
Acknowledgements
The authors thank the members of the New Brunswick Tree Improvement Council and Department of Natural Resources, Nova Scotia, for granting us access to their clonal orchards. The authors are grateful to D. Carleton and numerous summer students for data collection. M. Wright, K. Tosh, and M. Fullarton provided valuable orchard information. They also thank S. Heard, D. Gregoire, and D. Carleton for comments on earlier versions of this manuscript. Funding for this study was provided by the Atlantic Innovation Fund, Natural Science and Engineering Research Council of Canada, and the Growing Forward initiative.

