Killing two enemies with one stone? Genomics of resistance to two sympatric pathogens in Norway spruce

Trees must cope with the attack of multiple pathogens, often simultaneously during their long lifespan. Ironically, the genetic and molecular mechanisms controlling this process are poorly understood. The objective of this study was to compare the genetic component of resistance in Norway spruce to Heterobasidion annosum s.s. and its sympatric congener Heterobasidion parviporum. Heterobasidion root‐ and stem‐rot is a major disease of Norway spruce caused by members of the Heterobasidion annosum species complex. Resistance to both pathogens was measured using artificial inoculations in half‐sib families of Norway spruce trees originating from central to northern Europe. The genetic component of resistance was analysed using 63,760 genome‐wide exome‐capture sequenced SNPs and multitrait genome‐wide associations. No correlation was found for resistance to the two pathogens; however, associations were found between genomic variants and resistance traits with synergic or antagonist pleiotropic effects to both pathogens. Additionally, a latitudinal cline in resistance in the bark to H. annosum s.s. was found; trees from southern latitudes, with a later bud‐set and thicker stem diameter, allowed longer lesions, but this was not the case for H. parviporum. In summary, this study detects genomic variants with pleiotropic effects which explain multiple disease resistance from a genic level and could be useful for selection of resistant trees to both pathogens. Furthermore, it highlights the need for additional research to understand the evolution of resistance traits to multiple pathogens in trees.

. Although the understanding of major genes conferring disease resistance to single diseases in plants has advanced, the genetic and molecular mechanisms controlling quantitative disease resistance traits and its effectiveness against multiple attackers remains scarce, particularly in trees (Abdullah et al., 2017;Chen et al., 2018;Corwin & Kliebenstein, 2017;Ismael et al., 2020;Weiss et al., 2020).
Quantitative resistance traits have a continuous distribution of phenotypes from susceptible to resistant and are controlled by quantitative trait loci (QTL) -multiple loci with small to moderate effects (Corwin & Kliebenstein, 2017;Nelson et al., 2018). Quantitative disease resistance is assumed to be nonstrain specific and therefore durable (Ismael et al., 2020;Nelson et al., 2018;Wiesner-hanks & Nelson, 2016), however it is not always effective against different pathogens (Corwin & Kliebenstein, 2017). The nature of disease resistance to multiple pathogens could theoretically be explained from an organism level to a single gene level (Wiesner-hanks & Nelson, 2016). At the organism level, individuals can be resistant to multiple diseases because different unlinked QTLs present in an organism's genome are effective against different diseases independently (Risterucci et al., 2003;Wiesner-hanks & Nelson, 2016). At the genic level, multiple disease resistance could arise through the linkage of clusters of loci effective against single diseases (Schweizer & Stein, 2011) or by individual pleiotropic genes, where the same gene confers resistance to multiple diseases (Nelson et al., 2018;Wiesnerhanks & Nelson, 2016;Wisser et al., 2011).
The mapping and identification of QTLs is typically done through linkage mapping studies or genome-wide associations studies (GWAS) (Nelson et al., 2018). To guarantee the success of these experiments, they must be performed with high precision and comparable infection systems between pathogens, which is particularly challenging in forest systems (Ismael et al., 2020;Quesada et al., 2010). In recent years the knowledge of conifer genomics has improved vastly, which has allowed for more detailed studies on the genomic architecture of disease resistance traits Lind et al., 2014;Weiss et al., 2020). Within conifers, a well-studied pathosystem that allows for precise phenotyping is stem-and root-rot caused by members of the Heterobasidion annosum s.l species complex (Bodles et al., 2007;Chen et al., 2018;Dalman et al., 2013;Lind et al., 2014;Mukrimin et al., 2018;Skrøppa et al., 2015;Steffenrem et al., 2016).
Speciation in the Heterobasidion annosum s.l species complex began with a split between the ancestor of the pine-infecting species H. annosum s.s. and H. irregulare, and the ancestor of the nonpineinfecting species H. parviporum, H. abietinum, and H. occidentale Dalman et al., 2010). Species in the complex generally display sexual and somatic incompatibility and have different host ranges (Garbelotto & Gonthier, 2013). H. parviporum and H. annosum s.s., however, readily infect Norway spruce and share much of the Norway spruce distribution on the European continent ( Figure   S1; Chen et al., 2015;Dalman et al., 2010;Garbelotto & Gonthier, 2013;Niemela & Korhonen, 1998).
Norway spruce (Picea abies L. Karst) is a dominant conifer in boreal forests in Europe with a vast current population size . The sequencing of the Norway spruce genome and subsequent work has allowed the description of the species' evolutionary history and population structure Nystedt et al., 2013;Wang et al., 2020). Norway spruce is divided into three main domains, probably as a result of refugia through glaciation periods: a northern (Fennoscandian) domain ranging from Norway in the west to central Russia and two other domains in the Alps and Carpathians, with signs of main domain admixture-probably linked to recent expansion following the last glaciation period Li, 2020;Tsuda et al. 2016). Recent studies have described the genetics of wood properties, growth, phenology traits Milesi et al., 2019) and resistance to H. parviporum (Chen et al., 2018;Elfstrand et al. 2020).
Resistance to H. parviporum in Norway spruce is heritable (Chen et al., 2018;Lind et al., 2014;Steffenrem et al., 2016) and is characterized by many genes with relatively small effects on resistance . One QTL in PaLAR3, a gene involved in the synthesis of catechin and linked to H. parviporum resistance in Norway spruce, is known to respond to other stressors such as H. annosum s.s., the blue-stain fungus Endoconiophora polonica, and mechanical wounding (Danielsson et al., 2011;Hammerbacher et al., 2014;Nemesio-Gorriz et al., 2016). Therefore, we hypothesised that quantitative resistance to H. parviporum could provide multiple-disease resistance to other members of the H. annosum s.l. species complex.
In this study we measured disease resistance traits to H. annosum s.s. and H. parviporum in a well-characterized Norway spruce population part of the Swedish Norway spruce breeding programme Chen et al., 2018Chen et al., , 2019Lind et al., 2014;Milesi et al., 2019).
The programme is a result of phenotypic selection of trees across Europe based on growth, survival, stem quality and vitality, resulting in the inclusion of seven recognized Norway spruce genetic clusters in the current breeding population Haappanen et al. 2015;Milesi et al., 2019). We formulated the specific hypotheses that (i) Norway spruce has variation in its resistance traits to H. annosum s.s., (ii) resistance to H. annosum s.s. is correlated to resistance to H. parviporum, and (iii) QTLs could explain multiple-disease resistance in Norway spruce. To test these hypotheses, we studied resistance traits in 400 Norway spruce half-sib families following inoculation with H. annosum s.s. using quantitative genetics and genome-wide association methods (GWAS). Furthermore, we compared additive genetic resistance in half-sib families phenotyped for both H. annosum s.s and H. parviporum and identified potential multiple disease resistance QTLs with pleiotropic effects using multitrait GWAS.

| Plant material
A total of 400 open pollinated half-sib families from members of the founder population of the Swedish Norway spruce breeding programme were sown in 2016 (18 seedlings/family). After the first growth season, seedlings were randomised into a complete block design with three replications (Figure 1a), where each family was planted in 4-tree row-plots in plastic trays consisting of 24 separate 0.124 L plastic pots. The seedlings were grown for another season in Skogforsk's experimental nursery at Ekebo, Sweden (55°56′53.1″N 13°6′52.2″E) and subjected to standard watering and fertilisation.
No fungicides were used during cultivation.

| Artificial inoculations and phenotyping
Artificial inoculations were performed as described in Chen et al. (2018) with H. annosum s.s. Sä 16-4. The fungus was grown on Hagem's media plates for three weeks prior the experiment together with 5 mm P. abies wood plugs. Immediately prior to inoculation, bark was removed with a 6-mm diameter cork borer at 10 cm from the base of the seedling. A wooden plug colonised by H. annosum s.s. was then placed at the wound and covered with Parafilm (Chen et al., 2018). Ambient light and temperature conditions were maintained for 21 days, after which plants were harvested (from 20 August 2018 onwards).
Upon harvest, the diameter at the point of inoculation (D) was recorded and the lesion length (LL) above and below the edge of the inoculation point on the inner side of the bark was measured.
Sapwood growth of the fungus (SWG) was measured according to Arnerup and collaborators (2010): The inoculated stem was cut up into 5-mm discs and placed on moist filter paper in 9 cm Petri dishes together with the original colonised wooden plug. To avoid contamination, the stem was cut from the tip to, and the base to the point of inoculation, respectively. After seven days incubation under humid conditions, the presence of H. annosum s.s. on the discs was determined by observation of characteristic conidiophores under a stereomicroscope (Arnerup et al., 2010;Swedjemark et al., 1997).
Time of bud-set of seedlings following the first growing season, from mid-October to mid-November 2017, was recorded twice per week, with "1" and "0" representing the presence and absence of a visible bud, respectively.
Out of the 400 half-sib families phenotyped for H. annosum s.s., 269 were previously phenotyped for the same resistance traits to H. parviporum and reported by Chen et al. (2018).

| Statistical analyses
Measured traits were checked for recording errors and normality.
From a total of 5,924 observations, those with SWG =0 and no F I G U R E 1 Experiment set up. (a) Genotyping and phenotyping. Mother trees were genotyped and SNPs were filtered. Thereafter tree origin prediction and GWAS was performed. Half-sib families from the genotyped mother trees where phenotyped in three different blocks. These values were used to calculate EBVs, which in turn were used to calculate genetic correlations and the GWAS. (b) The half-sib families were phenotyped for resistance traits against H. annosum s.s. (N = 400) and H. parviporum (N = 501). The families phenotyped for both pathogens (N = 269) were used to calculate genetic correlations. Due to genotype filtering based in SNPs missingness, only a subset of mother trees met the cutoff and was used for origin assignment and GWAS conidiophores observed at either the point of inoculation or the inoculation plug were excluded from analyses (N = 235). Due to experimental errors progenies from the first block, with more than 75% of the seedlings scoring SWG =0, were also excluded (N = 69 observations). Resistance traits to H. parviporum phenotyped by Chen et al. (2018) were reanalysed in accordance with our criteria to remove bias. As LL showed a significant deviation from a normal distribution, the data was log-transformed, and a 0.5 constant was added to each value. Variance and covariance components for each trait were estimated using ASReml-R 4 (Butler et al., 2007) and the following linear mixed model was fitted for each trait individually: Where y ijkl is each observation on the lth seedling from the kth family in the jth block, is the general mean and B i is the fixed effect of the jth block. The variable F k is the random effect of the kth family, e ijkl is the random residual effect and D ijkl is a covariate for diameter at inoculation point. Wald tests were used to estimate the significance of fixed factors. Estimated breeding values (EBVs) for each family were defined as the coefficients of the random effect. Genetic correlation between traits was assessed by testing the association between EBVs using Pearson's product moment correlation in R.
The individual-tree narrow-sense heritability for each trait was estimated using the equation: where h 2 i ,̂ Where y ijklm is each observation on the lth seedling, at the mth week, from the kth family in the jth block where "1" corresponds to presence and "0" to absence of buds in the seedling, is the general mean and B i is the fixed effect of the jth block. The variable F km is the random effect of the repeated measurements for the kth family and G mk is the random effect of the mth week within the kth family, with a first order auto regressive variance assumption and e ijklm is the random residual effect. EBVs for each family were defined as the coefficients of the random effect of F k .

| SNP identification
Mother trees to the half-sib families were genotyped using 40,018 probes to cover intragenic regions in 26,219 P. abies gene models (Vidalis et al., 2018). DNA extraction, sequencing, and initial variant calling is described elsewhere Bernhardsson et al., 2020).

Variants
were filtered according to Bernhardsson et al. (2020) with minor modifications. Briefly, only biallelic
Missing variants in the remaining individuals were imputed with beagle 4.1 (Browning & Browning, 2007).

| Mother trees origin assignment
The ancestral origin of mother trees was assessed following Chen et al.
(2019) based on genotype similarity to individuals with known origin collected across P. abies natural range. Coordinates of the first five principle components of P. abies trees, from a sample population of 2572 (Li, 2020), with documented geographic origins and representative of the seven main genetic clusters were used as a training set in a "Random Forest" regression model ("randomForest" v4.6-14 package [Liaw & Wiener, 2002], r software v.3.3.1). The coordinates of the first five components of unknown individuals were then used to assign each mother tree to a given genetic cluster. The procedure was repeated 200 times with 8,000 iterations to estimate the accuracy of each assignment. Assignment of mother trees to a genetic cluster was determined to be true when the same allocation was repeated on more than 98% of occasions.

| SNP phenotype associations
Genome wide associations using different data sets were performed.
For H. annosum s.s., 330 mother trees were included after filtering for genotyping quality and relatedness (see above; Figure 1b). In order to perform multitrait GWAS between resistance traits to both H. annosum s.s. and H. parviporum we used the 220 overlapping mother trees between the population phenotyped for H. annosum s.s. resistance in this study and the population used in Chen et al. (2018) and Elfstrand et al. (2020). Associations were tested with GEMMA (Zhou & Stephens, 2012. EBVs calculated with ASreml R-4 (Butler et al., 2007) were used as the phenotype for each trait and kinship was accounted for with a standardized kinship matrix calculated in gemma (Zhou & Stephens, 2012. Principal component analysis (PCA) was computed with plink 1.9 (Chang et al., 2015) and used to identify and remove mother trees that were either too different or had very close family relationships with one another. Additionally, to account for population structure, three to four principal components were used as covariates depending on the subset of samples. Only Where y is a matrix of n × d traits, W a matrix of c × d covariates (fixed effects), is a matrix of the corresponding coefficients, x is an n -vector of the SNP genotypes, is a d vector of effet sizes for the d phenotypes, U is an n × d matrix of the random effects and is an n × d matrix of errors (Zhou & Stephens, 2012. Wald association tests were performed for each analysis testing the alternate hypothesis ≠ 0. In order to correct for multiple comparisons, False discovery rate (FDR) and Bonferroni, corrections were calculated with R. Since very few markers were significant following multiple comparisons correction, a suggestive significance threshold of 1x10 −5 (equivalent to the to 99.9 percentile) was used to identify candidate genes. The proportion of phenotypic variance explained by the SNP (PVE) was calculated according to (Shim et al., 2015).
The multitrait combinations were selected based on hypothesized relationships between traits, namely LL and SWG, within the experiment (different traits for the same pathogen) and between pathogens (same trait for different pathogens).   (Table 1).

| Resistance to H. annosum s.s. is not correlated to resistance to H. parviporum in Norway spruce
Narrow sense heritability estimates (h 2 ) were 0.49 for LL and 0.69 for SWG (Table 1) and a positive correlation between traits was observed ( Table 2) (Table 1). Correlation of the resistance traits in response to H. annosum s.s. and H. parviporum inoculations was low and nonsignificant (0.06 for LL and 0.08 for SWG; Table 2).
To test if there was a geographic effect on resistance, the ancestral origin of mother trees (i.e., before they were introduced in the Swedish breeding programme) was inferred based on genotype similarity to trees of known origin. One tree was assigned to the and northernmost clusters following a latitudinal cline (Figure 2), but that was not the case for SWG or any other phenotypes in H. parviporum ( Figure S2) (Table 3). Only eight markers could be placed in the linkage map and these were distributed in seven different linkage groups ( Figure S4). Interestingly, two of the SNPs detected specifically in Furthermore, one SNP (MA_99821:7939) was found within a gene annotated as an "ethylene responsive transcription factor" (Table 3).
A closer inspection of the gene model MA_99821g0010 shows that the gene indeed is a more likely orthologue of Cytokinin response factor 2 (AtCRF2) in Arabidopsis. Several SNPs in Pentatricopeptide repeat protein-and Tetraspanin genes were also detected (Table 3).
No QTLs were found to be associated with resistance to H. parviporum in this, or previous studies Mukrimin et al., 2018).

| Multitrait GWAS identifies loci with pleiotropic effects on resistance in Norway spruce
In order to test if loci have pleiotropic effects on the same trait for resistance to both pathogens, a multitrait GWAS was performed in 220 half-sib families. Considering the same significance threshold as above (p < 1.10 −5 ), 12 SNPs were found to be associated with LL and 7 with SWG ( Figure 3; Table 4). We then investigated correlations in allele effect sizes by plotting the effect sizes of all SNPs for resistance to H. parviporum as a function of their respective effect sizes in resistance to H. annosum s.s. (Figure 3). The SNPs were classified as belonging to two main categories (i) those with the same effect size positive pleiotropic effect and is co-located within 10 centimorgans (cM) from two different SNPs found to be significant in individual GWAS for SWG for both pathogens ( Figure S4;  (Table 2). This could be explained by differences in the pathogens' life strategy (Garbelotto & Gonthier, 2013;Hu et al., 2020;Oliva et al., 2011Oliva et al., , 2013, and by the ability of H. annosum s.s. to infect Pinus, using mechanisms that could also be effective when infecting Picea (Dalman et al., 2013) (Bartholomé et al., 2020). Consequently, the observed quantitative resistance to H. annosum s.s. and H. parviporum is likely to be dependent on both the environment in which infections take place and the genetic variation in resistance, which may have evolved independently to both Heterobasidion species.
The LL in response to H. annosum s.s. inoculation was significantly different in different genetic clusters of Norway spruce and followed a latitudinal cline; with mother trees from the Alpine domain having the longest lesions and trees from Southern and central Sweden being the most resistant in response to H. annosum s.s.
( Figure 1), but not to H. parviporum ( Figure S2). This is, to the best of our knowledge, the first time that a difference between tree origins has been observed in the interaction between a conifer and H. annosum s.l. (Bodles et al., 2007), although provenance effects on disease resistance have been reported for other forest pathogens (Hamilton et al., 2013;Perry et al., 2016). Moreira et al. (2014) observed that the level of constitutive defence in pines increases in species from higher latitudes and colder environments and is negatively correlated with early plant growth (Moreira et al., 2014). In Norway spruce quantitative traits such as growth and spring phenology follow environmental gradients in Europe   is worth noting that resistance traits in Norway spruce to H. parviporum are also correlated to the diameter at the inoculation point (Chen et al., 2018), but no significant difference between Norway spruce genetic clusters was observed in this interaction ( Figure S2). This is possibly influenced by the fungi´s respective tissue preferences, as H. annosum s.s. grows preferentially in the cambium and phloem tissues, while H. parviporum is concentrated in the sapwood and heartwood tissues (Hu et al., 2020;Oliva et al., 2011). An interaction located in the cambium and phloem tissues would be more susceptible to seasonal changes in fluxes, as shown previously in Norway spruce (Krokene et al., 2012).

| Novel gene models associated with resistance traits against H. annosum s.s.
Novel QTLs associated with resistance traits to H. annosum s.s. were found, four of which were exclusively found using multitrait associations ( Figure 3; Table 3). Recent use of multitrait GWAS in plant systems have proved useful in increasing the discovery power and understand the genetic make-up of complex traits such as response to stressors or leaf morphology (Chhetri et al., 2019;Thoen et al., 2017).
One advantage of this method is that the analysis of different traits jasmonate is the main hormonal pathway activated (Arnerup et al., , 2013Lundén et al., 2015), but recently the role of ABA has been highlighted (Kovalchuk et al., 2019). Because of the quantitative and potentially polygenic nature of the resistance traits in Norway spruce, it is likely that hormonal cross-talking takes place in the tissues in order to deploy a successful defence response.
Interestingly, other groups of SNPs in gene models associated nucleotide-binding site-leucine-rich repeat-type resistance genes (Källman et al., 2013). Here, we found a SNP in an argonaute1-like gene model associated to LL in both pathogens ( Figure 3). This gene model is the orthologue of Argonaute1 in Arabidopsis, which is known to modulate defence responses against bacterial and fungal pathogens by utilising endogenous small RNAs (Ellendorff et al., 2009;Katiyar-Agarwal et al., 2006). Interestingly this regulatory pathway is also utilised by pathogens like Botrytis cinerea, which use their own small RNAs via Argonaute1 to silence specific pathways in the host to establish successful infections (Weiberg et al., 2013). Given that pentatricopeptide repeat proteins and tetraspanins are also involved in RNA-mediated defence in Arabidopsis (Cai et al., 2018;Katiyar-Agarwal et al., 2006;Park et al., 2014) it is possible that candidate genes belonging to these groups, which were highlighted in this study, are involved in RNA mediated defence in Norway spruce against H. annosum s.s.

| Multitrait GWAS identifies pleiotropic QTLs associated with H. annosum s.s. and H. parviporum
Given that the resistance traits to H. annosum and H. parviporum had no correlation, it is not surprising that the SNPs associated with either pathogen in the univariate analysis were different (Table S1). It is worth mentioning that the exonic probes used cover only ~39% of the predicted gene models in the spruce genome (Vidalis et al., 2018) and that assembly of the genome is highly fragmented Nystedt et al., 2013). There could therefore, be significant variation associated to loci that are not observed in this study. Nonetheless, multitrait GWAS was used here to identify SNPs associated with resistance traits to both pathogens. A number of SNPs had effects that contributed to resistance traits to both H. annosum and H.
Interesting examples are the three different SNPs located within 10 cM in linkage group 3 ( Figure S4; Table S2). Two of these SNPs were found independently in the univariate GWAS for SWG for both pathogens and one other in the multitrait model for SWG  (Koutaniemi et al., 2015;Laitinen et al., 2017) which is associated with resistance to H. parviporum ). This gene TA B L E 4 SNPs associated to the same traits (lesion length (LL) and sapwood growth (SWG)) in both H.annosum s.s. and H. parviporum is specifically and differentially expressed in tissues after infection by H. parviporum, and is likely to be involved in the formation of the ligno-suberized boundary zone ).
Ligno-suberized boundary zone formation is a common feature of angiosperm and gymnosperm trees in response to a wide range of pathogens (Pearce, 1996;Woodward, 1992), which is in line with the synergistic pleiotropic effect observed in PaLAC5. Therefore, these results indicate that disease resistance to these two pathogens exists at genic level.
Another group of SNPs had the opposite effect for the same trait to the two pathogens (antagonist pleiotropy). Among the gene models harbouring such variants are an LRR-kinase receptor (MA_404302_2414) which is positively associated with resistance to H. annosum but negatively associated with resistance to H. parviporum (Figure 3, lower-right quadrant). LRR receptors with kinase functions are important components of both innate immunity and effector-triggered immunity in plants (Nürnberger & Kemmerling, 2006;Zhao et al., 2009). This particular LRR receptor harbours a conserved Malectin domain which is known to determine nonhost resistance in barley to powdery mildew strains adapted to wheat (Rajaraman et al., 2016). It is therefore possible that this LRR receptor recognises specific molecular patterns in only one of the pathogens leading to a successful defence response. Likewise, a secoisolariciresinol dehydrogenase-like gene (Figure 3, lower-right quadrant), which encodes for an enzyme involved in the production of matairesinol (Suzuki & Umezawa, 2007) had a negative pleiotropic effect. Lignans, such as matairesinol, have been shown to inhibit the activity of extracellular enzymes produced by a H. annosum s.l. isolate in vitro (Johansson et al., 1976;Popoff et al., 1975). In summary, SNPs associated to resistance traits to both pathogens can also have antagonistic pleiotropic effects on the infection outcome.

| Implications for disease resistance breeding in Norway spruce to Heterobasidion root-and stem-rot
Understanding the genetic architecture of tree resistance traits is an important task, as it will facilitate the development of resistance breeding strategies and ultimately ensure the success of reforestation programmes in the future (Buggs, 2020;Hall et al., 2016;Sniezko & Koch, 2017). H. annosum s.l. remains as one of the most devastating forest pathogens in the northern hemisphere and improved resistance to this species complex would be a desirable trait in breeding programmes (Garbelotto & Gonthier, 2013 show that some SNPs have a synergic pleiotropic effect, and selection based on these markers could be a useful strategy in breeding for resistance to both pathogens simultaneously. Furthermore, the significant variation in resistance to H. annosum s.s. with the predicted geographical origin of the mother trees indicates that disease resistance should be further studied in the ongoing assisted migration of Norway spruce trees.

| CON CLUS IONS
Here, we have used quantitative genetics together with exomecapture genomic data to understand the genetics behind resistance in Norway spruce to two closely related forest pathogens. The re- Furthermore, we found that these uncorrelated traits are associated with genomic variation in gene models with antagonist and synergic pleiotropic effects which are potentially involved in disease resistance, such as an PaLAC5, an LRR-kinase receptor and a secoisolariciresinol dehydrogenase. The QTLs with a synergic pleiotropic effect are an example of multiple disease resistance at the genic level and are of special interest as they could be utilised to select for trees with higher resistance to multiple pathogens. On the other hand, markers with an antagonistic pleiotropic effect could explain why these pathogens have evolved to inhabit different niches when infecting conifers. Finally, the results of this study highlight the need for further research to understand the plasticity of resistance traits in response to different pathogens under different environments -a key aspect in the success of reforestation programmes.

ACK N OWLED G EM ENTS
We would like to thank all the seasonal workers at the Skogforsk

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study have been