Genetic basis of growth, spring phenology, and susceptibility to biotic stressors in maritime pine

Abstract Forest ecosystems are increasingly challenged by extreme events, for example, drought, storms, pest attacks, and fungal pathogen outbreaks, causing severe ecological and economic losses. Understanding the genetic basis of adaptive traits in tree species is of key importance to preserve forest ecosystems, as genetic variation in a trait (i.e., heritability) determines its potential for human‐mediated or evolutionary change. Maritime pine (Pinus pinaster Aiton), a conifer widely distributed in southwestern Europe and northwestern Africa, grows under contrasted environmental conditions promoting local adaptation. Genetic variation at adaptive phenotypes, including height, spring phenology, and susceptibility to two fungal pathogens (Diplodia sapinea and Armillaria ostoyae) and an insect pest (Thaumetopoea pityocampa), was assessed in a range‐wide clonal common garden of maritime pine. Broad‐sense heritability was significant for height (0.219), spring phenology (0.165–0.310), and pathogen susceptibility (necrosis length caused by D. sapinea, 0.152; and by A. ostoyae, 0.021, measured on inoculated, excised branches under controlled conditions), but not for pine processionary moth incidence in the common garden. The correlations of trait variation among populations revealed contrasting trends for pathogen susceptibility to D. sapinea and A. ostoyae with respect to height. Taller trees showed longer necrosis length caused by D. sapinea while shorter trees were more affected by A. ostoyae. Moreover, maritime pine populations from areas with high summer temperatures and frequent droughts were less susceptible to D. sapinea but more susceptible to A. ostoyae. Finally, an association study using 4227 genome‐wide SNPs revealed several loci significantly associated with each trait (range of 3–26), including a possibly disease‐induced translation initiation factor, eIF‐5, associated with needle discoloration caused by D. sapinea. This study provides important insights to develop genetic conservation and breeding strategies integrating species responses to biotic stressors.


| INTRODUC TI ON
Forest ecosystems are challenged worldwide by changing environmental conditions (Turner, 2010). Warmer and drier climates are expected to increase fire risk, droughts, and insect outbreaks while warmer and wetter climates will probably increase storm and pathogen incidence on forests (Seidl et al., 2017), leading to episodes of high tree mortality (Allen et al., 2010) and severe economic losses (Hanewinkel et al., 2013). Changing environmental conditions can also cause range shifts in previously locally restricted pests and pathogens or shifts to increased pathogenicity (Desprez-Loustau et al., 2006). Thus, understanding forest tree genetic variation in disease (and other stressor) responses, in relationship to growth and spring phenology, is crucial to develop informed forest restoration, conservation, and management strategies. Moreover, genes underlying adaptive traits, such as growth or disease response, can serve tree breeding and increase forest productivity, for example, targeting resistance to drought or against pathogens in forest plantations (Neale & Kremer, 2011).
Forest trees are long-lived, sessile organisms. They are characterized by outcrossing mating systems, high standing genetic variation, large effective population sizes, and the production of vast numbers of seeds and seedlings exposed to strong selection (Petit et al., 2004;Petit & Hampe, 2006). High genetic and phenotypic differentiation has been observed in tree species along environmental gradients (e.g., Savolainen et al., 2007Savolainen et al., , 2013 or between contrasting habitats, indicating local adaptation (e.g., Lind et al., 2017). Common garden experiments (i.e., experiments evaluating trees from a wide range of populations under the same environmental conditions) provide valuable insights on the phenotypic and genotypic variation of forest trees (Morgenstern, 2011). They revealed genetic differentiation for adaptive traits (such as leaf phenology, i.e., flushing and senescence or growth) along latitudinal and altitudinal gradients (Mimura & Aitken, 2007;Vitasse et al., 2009). Geographical variation has also been found for disease resistance against certain pests (Menéndez-Gutiérrez et al., 2017) and pathogens (e.g., Freeman et al., 2019;Hamilton et al., 2013). Phenological traits, such as time of flowering, or spring bud burst and autumn leaf senescence, are often genetically correlated with disease susceptibility, providing hints on the different kind of resistance or avoidance mechanisms found in forest trees (Elzinga et al., 2007).
Disease resistance is generally thought to be the result of selective pressures exerted by the pathogen, in areas where host and pathogen have co-existed during considerable periods of time, under the co-evolution hypothesis (e.g., Burdon & Thrall, 2000;Ennos, 2015). In this line, geographical variation in disease resistance has been interpreted in some cases as a result of past heterogeneous pathogen pressures within the range of a given host species (Ennos, 2015;Perry et al., 2016). However, the past distribution of pathogen species is often unknown (Desprez-Loustau et al., 2016); therefore, other processes than co-evolution, such as "exaptation" (Gould & Vrba, 1982) or "ecological fitting" (Janzen, 1985), should not be excluded. These biological processes have been suggested when, for example, variability in disease resistance was observed in tree species with no co-evolutionary history with a pathogen (Freeman et al., 2019;Leimu & Koricheva, 2006;Newcombe, 1996). Exaptation refers to an adaptive trait, for example, resistance to a novel, possibly invasive pathogen, that evolved in response to another selective pressure. Similar to this, in the case of "ecological fitting," tree species that expanded their range encounter previously unknown pathogens in the new environment (Agosta & Klemens, 2008). Disease resistance in these cases may have evolved in response to other pathogens but had broad-range efficacy, even to a novel pathogen.
Generic mechanisms of resistance in conifers include the production of large amounts of non-volatile compounds (resin acids) that can act as mechanical barriers to infections (Phillips & Croteau, 1999;Shain, 1967) and volatile compounds (such as monoterpenes or phenols) that can be toxic to fungi (Chou & Zabkiewicz, 1976;Cobb et al., 1968). The composition of secondary metabolites can show marked geographic variation in tree species (Meijón et al., 2016). The evolution of plant defenses against biotic stressors can also be shaped by differences in resource availability and environmental constraints throughout the host's species distribution. Depending on resource availability, plants have evolved distinct strategies by investing either more in growth, to increase competition ability, or more in chemical and structural defenses, to better respond to herbivores and pathogens (Herms & Mattson, 2004). Typically, faster-growing trees invest more in inducible defenses while slow-growing trees invest more in constitutive defenses (Moreira et al., 2014).
Many quantitative traits in forest species, including height, spring phenology, and disease resistance, show significant heritability and often stronger differentiation (Q ST ) among populations than neutral genetic markers (F ST ) (Hamilton et al., 2013; see review in Lind et al., 2018), a common indication of adaptive divergence. Major-effect genes for growth, for example, korrigan in Pinus pinaster (Cabezas et al., 2015), and resistance genes against forest pathogens, for example, against the fusiform rust disease in Pinus taeda (Kuhlman et al., 2002) and against white pine blister rust in several other North American pine species (Sniezko et al., 2008;Weiss et al., 2020), have been identified in forest trees. However, most adaptive traits have a highly polygenic basis of quantitative inheritance, typically involving many loci with rather small effects (De la Torre et al., 2019;Goldfarb et al., 2013). Most association genetic studies in forest trees focused on wood property and growth traits to assist tree breeding (e.g., Beaulieu et al., 2011;Neale et al., 2006;Pot et al., 2005). In addition, genetic association approaches identified some loci associated with K E Y W O R D S association genetics, genetic correlations, heritability, pathogen susceptibility, Pinus pinaster, spring phenology other ecologically important traits, such as cold hardiness (e.g., Eckert et al., 2009;Holliday et al., 2010), drought tolerance (reviewed in Moran et al., 2017), or disease resistance (e.g., Liu et al., 2014;Resende et al., 2017). Compared with more intensively studied traits, association studies addressing biotic interaction traits, including responses to pests and pathogenic fungi, are still limited to few tree species.
Our study focused on maritime pine (Pinus pinaster Aiton, Pinaceae), a long-lived conifer from southwestern Europe and northern Africa with a wide ecological amplitude. We assessed height, spring phenology (bud burst and duration of bud burst as phenological phases of tree growth), and susceptibility to pests/ pathogens in a clonal common garden, which allowed us to explore variation in disease response and genetic correlations with other traits in range-wide populations of maritime pine. Considering disease and growth traits together is relevant from an evolutionary and ecological perspective and can have important implications in terms of management, especially in breeding programs. We selected three important disease agents: two fungal pathogens, Diplodia sapinea (Botryosphaeriaceae) and Armillaria ostoyae (Physalacriaceae), as well as the pine processionary moth, Thaumetopoea pityocampa (Thaumetopoeidae), a main defoliator of pine forests.
Maritime pine has a highly fragmented natural range in the western Mediterranean Basin, the Atlas Mountains in Morocco, the Atlantic coast of southern France and the west coast of the Iberian Peninsula, and grows from sea level to over 2000 m altitude. Genetic diversity and population structure are high in natural populations of maritime pine, especially in the Iberian Peninsula, possibly due to its long-term persistence in this region (Bucci et al., 2007;Jaramillo-Correa et al., 2015;Petit et al., 1995;Salvador et al., 2000). In addition, traits, such as stem form, height (González-Martínez et al., 2002), metabolite content (Meijón et al., 2016), drought response (Aranda et al., 2010;Gaspar et al., 2013), and pest and disease resistance (Burban et al., 1999;Desprez-Loustau & Baradat, 1991;Elvira-Recuenco et al., 2014;Schvester, 1982), are highly variable in maritime pine and often strongly differentiated across geographical regions. Maritime pine has also been widely planted and is currently exploited for timber and paper, for example, covering ~1.03 million ha in southwestern France (including the Landes region), one of the largest plantation forests in Europe (Memento FCBA 2020, https://www.fcba.fr/resso urces/ memen to-2020/). Despite the ecological and economic importance of maritime pine natural forests and plantations, only a few genetic association studies have been developed in this species. Lepoittevin et al. (2012) identified two loci associated with growth and wood cellulose content, respectively. Cabezas et al. (2015) revealed four SNPs in korrigan (a pine gene orthologous to an Arabidopsis degrading enzyme cellulase) also significantly associated with growth traits (total height and polycyclism). Finally, Bartholomé et al. (2016) reported four loci for stem straightness and three loci for height growth. In addition, Budde et al. (2014) were able to predict 29% of the phenotypic variation in a fire adaptive trait (proportion of serotinous cones) in eastern Spain based on 17 SNP loci. However, none of these studies targeted spring phenology or biotic interaction traits, such as disease resistance.
Diplodia sapinea is the causal agent of several diseases, such as tip-blight, canker or root collar necrosis in needles, shoots, stems, and roots of conifers, eventually leading to mortality in case of severe attacks (Luchi et al., 2014;Piou et al., 1991). The pathogenicity of D. sapinea is associated with environmental conditions. It can remain in an endophytic form, that is, without causing any symptoms, until stressful environmental conditions, such as drought (Desprez-Loustau et al., 2006;Stanosz et al., 2002), hail storms (Zwolinski et al., 1990), or changes in the nitrogen concentration of the soil (Piou et al., 1991;Stanosz et al., 2004), weaken the host and trigger D. sapinea pathogenicity. Trees from all ages are affected (Chou, 1978;Georgieva & Hlebarska, 2017), though seedlings and old trees show increased susceptibility (Swart & Wingfield, 1991). The fungus can be found in many conifers, especially in the genus Pinus. Iturritxa et al. (2013) classified maritime pine as moderately susceptible in a comparative study of six pine species. The fungus was first described in Europe in 1823 under the name Sphaeria sapinea and then received many synonyms (Piou et al., 1991). Recent surveys showed that D. sapinea is currently broadly distributed in all pine forests throughout the world, though its origin is debated (Adamson et al., 2021;Brodde et al., 2019;Burgess et al., 2004). Serious damage associated with D. sapinea in Europe has only been reported in the last decades but it may become a serious threat to pine forests, as climate change will certainly favor pathogen activity by increasing temperature, and the frequency and intensity of drought events (Boutte, 2018;Desprez-Loustau et al., 2006;Woolhouse et al., 2005). In this line, recent outbreaks associated with D. sapinea in northern Europe suggest an ongoing northward expansion (Brodde et al., 2019).
Armillaria ostoyae is a root pathogen that causes white rot and butt rot disease in conifers, leading to growth deprivation, high mortality, and major losses in timber wood, hence its economic importance (Cruickshank, 2011;Heinzelmann et al., 2019). The species can be traced back to six million years ago, both in Eurasia and in North America (Koch et al., 2017;Tsykun et al., 2013). Armillaria ostoyae has been reported in all the conifer forests of the Northern Hemisphere but it appears to be replaced by A. mellea (Marxmüller & Guillaumin, 2005) in the Mediterranean due to higher temperatures and drought. It is likely that A. ostoyae would have co-existed for a long time with maritime pine in Europe (Tsykun et al., 2013), and consequently, it could have been affected by the same extinctionrecolonization events associated with past climatic changes as its host . Armillaria ostoyae is one of the most common fungal species in maritime pine forests, being particularly dangerous as it can act as both a parasite and a saprophyte (Cruickshank et al., 1997;Labbé, Lung-Escarmant, et al., 2017), that is, the death of its host does not prevent its spread. In maritime pine, the severity of A. ostoyae symptoms is related to host age, with higher mortality in young trees (Labbé et al., 2015;Lung-Escarmant & Guyon, 2004). Climate change is predicted to foster an increased impact of A. ostoyae on conifer forests in the coming years (Kubiak et al., 2017).
Thaumetopoea pityocampa is considered the most severe defoliator insect in pine forests in southern Europe and northern Africa (Jactel et al., 2015), leading to severe growth loss (Jacquet et al., 2013). The species typically reproduces in summer followed by larval development during autumn and winter. Caterpillars and moths of T. pityocampa are sensitive to climatic and environmental conditions, and this pest is expected to expand its range following events of climate warming (Battisti et al., 2006;Toïgo et al., 2017).
The specific objectives of our study are to (1) estimate phenotypic variability and heritability within and among range-wide populations of maritime pine for height, spring phenology, and susceptibility to pests/pathogens; (2) test for adaptive divergence across the maritime pine range for these traits (i.e., Q ST vs. F ST approach); (3) analyze the pattern of trait correlation within and among populations, in particular to identify putative correlations that could be useful/detrimental for conservation and breeding programs; and (4) identify loci associated with disease-related, height, and spring phenology traits using genotype-phenotype association. Altogether, our approach, combining the evaluation of a clonal common garden and a genotyping array, produced relevant insights on the evolution, genetic basis, and architecture of adaptive traits in maritime pine, an ecologically and economically important forest tree species.   Jaramillo-Correa et al., 2015). Open-pollinated seeds were collected from approximately 30-40 trees in each stand, separated at least 50 m from each other.

| Plant material and common garden measurements
One single seed from each tree was planted in a nursery. A selection of the saplings (the "mother plants") was propagated by cuttings to obtain genetically identical replicates. These clones (average of over 12 distinct genotypes per population) were used to establish the clonal common garden. The common garden design consisted of 517 clones (i.e., genotypes) planted in eight randomized complete blocks, with one clonal copy (ramet) of each genotype replicated in each block (total number of trees 4136). For the pathogen inoculation experiments, we chose genotypes from populations representing each of the six gene pools. Because of the higher logistical effort involved in inoculation protocols with respect to other measurements (see below), it was not possible to include all common garden genotypes in these experiments (see Table S1.1).
Height, spring phenology (including two growth phenophases, namely bud burst and duration of bud burst), and incidence of processionary moth (T. pityocampa) were measured in all individuals from 5-8 blocks, depending on the trait (sample size of 1,436-3,310 trees, see Table S1.1). Pathogen susceptibility was assessed in a subset of genotypes using excised branches collected from the clonal trial (sample size of 174-453 branches, see Table S1.1 and below). begin to space, 4: appearance of the needles, 5: elongation of the needles; see Figure S2.1). The Julian day of entry in each stage (S1 to S5) was scored for each tree. Julian days were converted into accumulated degree-days (0°C basis) from the first day of the year, to take into account the between-year variability in temperature. The number of degree-days between stages 1 and 4 defines the duration of bud burst. The presence or absence of pine processionary moth nests (T. pityocampa) in the tree crowns was assessed in March 2018.

| Experimental evaluation of necrosis length and needle discoloration caused by Diplodia sapinea
Direct inoculations on the trees in the clonal common garden were not possible, as this approach would be destructive and not compatible with the objective of long-term monitoring of the trees in the experimental garden. Therefore, we used detached branch inoculations as proxies to assess D. sapinea susceptibility of genotypes evaluated as necrosis length and needle discoloration. Inoculations were carried out in the laboratory on excised shoots taken from pines in the common garden (for a detailed laboratory protocol, see Supporting Information S3.1). We used the pathogen strain Pier4, isolated from P. nigra cones in Pierroton, France (May 2017), and maintained on malt agar medium. The identity of this strain as D. sapinea was confirmed by sequencing the ITS region, amplified using the primers ITS1-F and ITS4 (Gardes & Bruns, 1993), and blasting it against the NCBI nucleotide database (Benoît Laurent, personal communication). Only current-year shoots at the phenological stages 3 to 5, that is, with fully elongated buds but not fully mature, were collected, otherwise sampled randomly (see Supporting Information S2.1). We ensured that no disease symptoms, especially non indicating any attack by D. sapinea, were observed on the sampled trees. For the inoculation, we removed a needle fascicle in the middle of each shoot with a scalpel. A 5-mm-diameter plug of malt agar taken at the active margin of a D. sapinea culture was placed on the wound, myceliumside down, and then wrapped in cellophane. A total of 24 randomly chosen control shoots, one for each inoculation day, were treated in the same manner but with plugs of sterile rather than colonized malt agar. The shoots were put in water and kept in a climatic chamber set at 20°C with a daily cycle of 12 h of light and 12 h of dark (Blodgett & Bonello, 2003;Iturritxa et al., 2013). Six days after the inoculation, we removed the cellophane and measured the lesion length around the inoculation point with a caliper. The shoots were not lignified and the lesions were visible. However, the surface was superficially stripped to see the limit of the lesion when it was not visible otherwise. Needle discoloration was also observed and evaluated using a scale from 0: no discoloration to 3: all needles along the necrosis showed discoloration (see Figure S3.1). To confirm that discoloration was caused by the pathogen, one discolored needle from one branch per population was placed on a malt agar Petri dish to grow.
After three days, D. sapinea could be visually identified in each Petri dish. Although we could not rule out that D. sapinea might have been present as endophyte in the trees at the beginning of the experiment, branches treated with sterile malt agar plugs as control did not show any necrosis nor needle discoloration. This confirmed that our treatment did not trigger any pathogenic outbreak independent from the mycelium inoculations.
For this experiment, we sampled a total of 453 branches, from 151 genotypes (i.e., one branch from each of three replicate trees per genotype) in ten populations, representing all differentiated gene pools known in maritime pine (see Jaramillo-Correa et al., 2015). On 24 days, between June 12 and July 31, 2018, we collected between 10 and 40 branches, depending on technical schedule and availability. One lateral branch per tree was cut from the previous year whorl of the selected trees (see above) and taken to the laboratory for inoculation. Inoculations were performed on the leader shoot of the current-year whorl of the excised branch.

| Experimental evaluation of necrosis length caused by Armillaria ostoyae
Root inoculations in the clonal common garden were not possible.
Therefore, we used detached branch inoculations as proxies to evaluate A. ostoyae susceptibility measured as necrosis length. As common garden trials of forest tree species are very valuable and have been established and maintained over many years, the use of detached branches in experimental settings as proxies is a common practice to avoid destructive sampling. Using excised branches was shown to be a reliable approach, for example, for measuring cavitation resistance (e.g., Lamy et al., 2014), photosynthesis parameters and stomata conductance (Akalusi et al., 2021) and characterizing the timing of bud burst under experimental conditions (e.g., Viherä-Aarnio et al., 2014), as well as for evaluating the susceptibility to root pathogens (Bodles et al., 2007;Matusick et al., 2016;Nagy et al., 2006;Santos et al., 2015;Shain, 1967;Swedjemark & Karlsson, 2006;Underwood & Pearce, 1992), similarly to our study. Usually, excised branches maintain their "physiological attributes" after cutting at least during some hours or even days, especially in conifers most likely due to their xeromorphic tissue (e.g., Richardson & Berlyn, 2002;Warren, 2006). Although this approach is not ideal, it was the only way to obtain phenotypes related to infection with A. ostoyae on the trees from the clonal common garden. Alternatively, seedlings or potted saplings could have been used to evaluate pathogen resistance but due to ontogenic effects such an approach would also not be ideal.
For the inoculation with A. ostoyae, we used the pathogen strain A4, collected from a dying maritime pine tree in La Teste (Gironde, France) in 2010 (Labbé, Lung-Escarmant, et al., 2017).
For the experiment, two plugs of 5 mm diameter of malt agar with the A. ostoyae mycelium were placed on the top of a mixture of industrial vegetable soup (Knorr 9 légumes©, Heilbronn, Germany), malted water, and hazelnut wood chips in a 180-mL plastic jar (Heinzelmann & Rigling, 2016) (for a detailed laboratory protocol, see Supporting Information S3.2). The lid was closed loosely enough to allow some oxygen flow. The jars were placed during three months in the dark, in a heat chamber set at 23°C and 80% humidity before inoculation. The basal part of the sampled shoots (ca. 8 cm) was placed in the center of the mycelial culture in the heat chamber, maintaining the same temperature and humidity settings as for the mycelium growth, but adding an additional 12-h cycle of light/dark. Only the jars showing a minimum jar occupation by A. ostoyae of 60% were used. After three weeks, wood colonization success for each sample was evaluated visually by confirming the presence of mycelium under the bark. At this point, the needles and the phloem of the branches were still green and thus the sample was considered as physiologically alive and adequate for susceptibility phenotypic measurements. The length of the colonizing mycelium and length of the lesion in the sapwood (i.e., wood browning, hereafter referred to as necrosis) were measured. In the jar, we also visually evaluated the level of humidity of the medium (dry, medium, and very humid) and A. ostoyae growth.
A total of six control branches, one genotype randomly chosen from each of the six gene pools, were prepared in the exact same manner, but with plugs of sterile malt agar as opposed to those colonized by A. ostoyae.
For this experiment, we randomly sampled ten (except for Tamrabta with eight) maritime pine genotypes from six populations, one for each of the six differentiated gene pools represented in the CLONAPIN common garden. We selected fully elongated current-year shoots (bud burst stage 4 and 5), with a maximum diameter of 15 mm and a minimum length of 10 cm, as proxies.
A total of 174 branches from 58 genotypes (i.e., one branch from each of three replicate trees per genotype) were measured, cut, and taken to the laboratory to be inoculated, on October 3-4, 2018.

| Climatic data
Summary climate data for the years 1950-2000 were retrieved for 32 variables from Worldclim (Hijmans et al., 2005) and a regional climatic model (Gonzalo, 2007) (Eveno et al., 2008;Grivet et al., 2011), significant environmental associations with climate at the range-wide spatial scale (Jaramillo-Correa et al., 2015), or differential expression under biotic and abiotic stress in maritime pine . After standard filtering followed by removal of SNPs with uncertain clustering patterns (visual inspection using GenomeStudio v. 2.0), we kept 5176 polymorphic SNPs, including 4227 SNPs with a minor allele frequency (MAF) above 0.1.

| Quantitative genetic analyses
To estimate the genetic variance components of the analyzed traits, we ran two different sets of mixed-effect models. The first one included three hierarchical levels, to account for the strong population genetic structure in maritime pine: gene pool, population nested within gene pool, and genotype (clone) nested within population and gene pool, as follows: where for any trait y igjk , µ denotes the overall phenotypic mean, block i represents the fixed effect of experimental block i, gp g is the random effect of gene pool g, gp(pop) gj is the random effect of pop j nested within gene pool g, gp(pop(genotype)) gjk is the random effect of genotype k nested within population j and gene pool g, and ε is the overall residual effect. In model 2, cov represents the covariates implemented when modeling the presence of pine processionary moth nests (i.e., tree height in 2015) and necrosis caused by A. ostoyae (i.e., a categorical evaluation of jar humidity). Using these models, the effects at the gene pool level had extremely wide confidence intervals, probably because some of them were represented by only a few highly contrasted populations. Thus, we also ran a set of models without the gene pool effect, as follows: This set of models produced more accurate BLUP estimates, avoiding also confounded effects between gene pool and population, and was thus preferred for some of the subsequent analyses (see below).
All mixed-effect models were fitted in a Bayesian framework processionary moth nests and needle discoloration caused by D. sapinea infection that followed a binomial distribution and were modeled with logit and probit link functions, respectively. Multivariate-normal prior distribution with mean centered around zero and large variance matrix (10 8 ) was used for fixed effects with the exception of the model for needle discoloration caused by D. sapinea where a gelman prior for V was set, as suggested by Gelman et al. (2008) for ordinal traits. Inverse Wishart non-informative priors were used for the variances and covariances with a matrix parameter V set to 1 and a parameter n set to 0.002 (Hadfield, 2010). Parameter expanded priors were used to improve the convergence and mixing properties of the chain as suggested by Gelman (2006) for models on the presence of pine processionary moth nests, needle discoloration caused by D. sapinea, and necrosis caused by A. ostoyae. Parameter estimates were not sensitive to change in the priors. The models were run for at least 750,000 iterations, including a burn-in of 50,000 iterations and a thinning interval of 500 iterations. Four chains per model were run to test for parameter estimates convergence. Gelman-Rubin criterion potential scale reduction factor (psrf) was consistently below 1.01 (Gelman & Rubin, 1992) (see Table S4.1, for further details on model specifications).
Variance components from the first set of models (including gene pool effect) were then used to compute broad-sense heritability (H 2 ) as: where 2 genotype is the variance among genotypes within populations and gene pools and 2 e the residual variance. In addition, we computed the variance ratios associated with the population (H 2 p ) and gene pool (H 2 gp ) effects as: where 2 pop is the variance among populations within gene pools, and 2 gp is the variance among gene pools. When appropriate, we included an extra term in the denominator to account for implicit logit and probit link function variance (π 2 /3 and +1, respectively; Nakagawa & Schielzeth, 2010). Genetic differentiation among populations for the analyzed traits (Q ST ) was calculated as in Spitze (1993), using the second set of models (without gene pool effects, equations 3 and 4), as we were interested in the differentiation at the population level: where 2 pop is the overall variance among populations, and 2 genotype is the variance among genotypes within populations. Finally, we estimated the global F ST using all available SNP genotypes in SPAGeDi 1.5 (Hardy & Vekemans, 2002). The difference between global F ST and Q ST values for each adaptive trait was considered significant when the 95% credible intervals (CI) did not overlap.
Genetic correlations were calculated with Pearson's coefficient of correlation using the best linear unbiased predictors (BLUPs) of the genotype (clone) effect (Henderson, 1973;Robinson, 1991) for each trait. Genetic correlations reflect the proportion of variance shared by two traits and indicate pleiotropic effects of genes or linkage between loci affecting both traits. Genetic correlations are based on individuals from the same population, that is, which have experienced the same evolutionary history, and are relevant to identify trade-offs between traits (e.g., negative genetic correla-

| Genotype-phenotype association
First, we used a mixed linear regression approach (MLM, Yu et al., 2006) implemented in Tassel v. 5.0 (Bradbury et al., 2007) to identify single SNPs associated with each of the phenotypes. Phenotypic values were estimated using the BLUPs accounting for both population and genotype effects from the models without gene pool effects (second set, equations 3 and 4). Ancestry proportions of each sample to six genetic clusters (see Jaramillo-Correa et al., 2015) were computed using STRUCTURE (Pritchard et al., 2000). These ancestry proportions were included as covariates in the MLM. A covariance matrix accounting for relatedness between all sample pairs was estimated using Loiselle's kinship coefficient (Loiselle et al., 1995) in SPAGeDi 1.5 and was included as random effect. Negative kinship values were set to zero following Yu et al. (2006). Only loci with a pvalue below 0.005 in the Tassel analyses and with a minimum allele frequency above 0.1 were used for further analyses. Second, for the identified associations, we used a Bayesian mixed-effect association  Budde et al., 2014). As in the previous models, the STRUCTURE ancestry proportions were used as covariates and the relatedness matrix as random factor. Mean allelic effects (γ) and 95% confidence intervals were obtained from the distribution of the last 20,000 iterations (50,000 in total). Only those SNPs with credible intervals not overlapping zero were considered to have a significant (non-zero) effect on the trait.
Functional annotations, SNP motives, and BLAST results for significant genotype-phenotype associations were retrieved from Plomion et al. (2016). Finally, in order to inspect visually geographical patterns, the minimum allele frequency of significantly associated SNPs was estimated in each population using SPAGeDi 1.5 and plotted in a map.  Most traits showed strong differences among gene pools, although the variance ratio associated with the gene pool effect was characterized by wide confidence intervals (

| Correlations between traits and with environmental variables
All genetic correlations (i.e., those considering only the genotype effect) involving height and spring phenology were significant while Abbreviations: bb, bud burst; dbb, duration of bud burst; disc, needle discoloration; processionary, presence/absence of processionary moth nests; dd, degree-days.
Heritability and variance ratios for incidence of the processionary moth were computed using height as a covariate. For A. ostoyae necrosis length, the population and gene pool levels were identical as only one population per gene pool has been analyzed. Values in bold are significant. Values in squared brackets indicate the 95% credible intervals. disease susceptibility traits from excised branches did not show significant genetic correlations except for a correlation between necrosis length and needle discoloration caused by D. sapinea ( Table 2).
The strongest genetic correlation was observed between bud burst and duration of bud burst in 2015 (0.798, p-value < 0.001). Once the population effect was considered, a significant negative trait covariation was found between necrosis length caused by each of the two fungal pathogens (−0.692, p-value < 0.001; Table 2, Figure 3).
We also observed significant across population correlations with height, negative for necrosis length caused by A. ostoyae (−0.653, pvalue < 0.001), and positive for necrosis length caused by D. sapinea (0.679, p-value < 0.001).
Moreover, using mixed linear models that account for gene pool as random factor, we observed a negative relationship for maximum temperatures in the summer months with height, duration of bud burst in 2017, and necrosis length and needle discoloration caused by D. sapinea (Table 3)

| Genotype-phenotype association
Between three and 28 SNPs were significantly associated with each of the phenotypic traits evaluated under different genotype effect models (see Table S6.1). Here, we only focus on SNPs that were significant under the additive genetic model, this model being built on three genotypic classes and therefore considered the most robust.
Based on this model, five SNPs were associated with height, 27 SNPs were associated with spring phenology (considering altogether the different phenology traits and measurement years), and seven with pathogen susceptibility (Table 4). In total, four significantly associated SNPs showed non-synonymous changes. Two non-synonymous SNPs were associated with bud burst in 2017 ( Figure S6.1). In addition, one non-synonymous SNP was associated with needle discoloration caused by D. sapinea and another one to duration of bud burst in 2015 (Table 5 and Figure 5). All the remaining SNPs involved in significant associations under the additive model were either noncoding or the effect of the substitution was unknown (Table S6.1).
The geographical distribution of minor allele frequency of the associated SNPs was quite variable and did not reflect the population genetic structure of the species ( Figure S6.1-6.4).

| DISCUSS ION
In the current context of climate change, understanding the genetic basis of adaptive traits in tree species is key for an in-

| Genetic variation and differentiation, and correlations with climate
Broad-sense heritability estimates were in the range of previously published values in maritime pine or other forest tree species.
The variance ratios associated with population and gene pool effects reflect genetic differences that are not necessarily due to selection but might also reflect other evolutionary processes at   Note: Bayesian mean SNP effects and 95% credible intervals (CIs) were obtained from the distribution of the last 20,000 iterations in BAMD. Marker codes and linkage groups as reported in Plomion et al. (2016). Abbreviations: bb, bud burst; dbb, duration of bud burst; disc, needle discoloration; LG, linkage group; MAF, minimum allele frequency; unk, unknown; nc, noncoding; non-syn, non-synonymous; syn, synonymous.

TA B L E 4 (Continued)
different evolutionary histories, our results suggest that selection and hence adaptation at quantitative traits take place at a finer spatial scale. Therefore, genetic differentiation for adaptive traits at the population level (e.g., as estimated by Q ST ) is expected to be more informative. Neutral genetic differentiation, that is, F ST , was moderate (F ST = 0.109 [0.013; 0.325], p-value < 0.001) and significantly lower than Q ST estimates obtained for height and necrosis length caused by D. sapinea, which suggested that divergent selection is promoting local adaptation in these traits (Lamy et al., 2011;Whitlock & Guillaume, 2009).
Height is a crucial, frequently studied trait in forest trees, as increased height growth is a main target of breeding programs (e.g., Cornelius, 1994;Kremer & Lascoux, 1988). In our study, we found   ing the variation across populations. This is not surprising, as a big part of the trait variation was due to population effects, reflecting differences in their evolutionary history due to distinct levels of genetic drift and gene flow, as well as selective pressure (Archambeau et al., 2020). Height is known to be a highly integrative trait closely responding to abiotic factors, such as climate Jaramillo-Correa et al., 2015), and has thus been used in combination with genetic markers to identify relevant conservation units in maritime pine . In our study, we showed that height is not only correlated with climate but also with biotic factors, such as pathogen susceptibility indicated by necrosis length on excised branches (positively for D. sapinea and negatively for A. ostoyae). This is crucial information, for example, for the Landes maritime pine breeding program, which is based on artificial selection for height growth and straightness, and that may have selected at the same time for higher susceptibility to D. sapinea, a pathogen that is expected to increase its presence in maritime pine productive forests in the near future (see below).

Spring phenology traits (bud burst and duration of bud burst)
showed low-to-moderate broad-sense heritability, depending on the year ( (Ducousso et al., 1996), and are probably related to the higher risk of late frost in northern populations as it is the case in the French Atlantic range of maritime pine (see also, e.g., Alberto et al., 2013;Zohner & Renner, 2014). Such clines in spring phenology along climate gradients suggested local adaptation due to divergent selection, which might suffer mismatch with changing climate conditions (Alberto et al., 2013;Badeck et al., 2004;Lindner et al., 2010). Spring phenology can also play a role in resistance to or avoidance of forest tree pathogens (e.g., Ghelardini & Santini, 2009;Nielsen et al., 2017;Swedjemark et al., 1998). In line with this, we found a positive correlation across populations between needle discoloration and necrosis length caused by D. sapinea with spring phenology indicating that earlier flushing trees with faster spring growth showed less severe disease symptoms. Krokene et al. (2012) showed that the concentrations of starch and total sugars (glucose, fructose and sucrose) in twigs of Picea abies change during shoot development, which affects pathogen-related symptoms. In our study, inoculations were carried out on twigs with already elongated needles; however, the chemical composition of twigs might differ with time since bud burst and we cannot fully discard these effects in our assessment of disease incidence.  (Krokene et al., 2012). Additionally, detached branches were used in order to avoid destructive sampling in the clonal common garden. As mentioned in material and methods, pathogen inoculations on detached branches are commonly used, but can only represent a proxy to evaluate pathogen susceptibility, especially for root pathogens (Underwood & Pearce, 1992 et al., 2021) or pathogen pressure was not strong enough, differences in susceptibility among maritime pine populations might be due to exaptation, that is due to correlated traits selected for other functions (Agosta & Klemens, 2008). Populations of maritime pine strongly vary geographically in many traits related to growth and response to drought, along the gradient from North Africa to the Atlantic regions of the Iberian Peninsula and France (Aranda et al., 2010;Corcuera et al., 2012;Correia et al., 2008;Gaspar et al., 2013;de la Mata et al., 2014). Some of these traits may indirectly influence their susceptibility to pathogens, as observed here for D. sapinea. For example, faster-growing maritime pine trees from northern populations are known to invest more in inducible defenses while slow-growing trees from southern populations invest more in constitutive defenses (López-Goldar et al., 2018). The positive correlation between height and necrosis length caused by D. sapinea might indicate that constitutive defenses confer better resistance to this pathogen in the southern populations.
Also, Meijón et al. (2016) (Labbé, Lung-Escarmant, et al., 2017;Prospero et al., 2004;Solla et al., 2011), as the pathogen grows faster once it enters the organism and reaches the cambium (Solla et al., 2002). Infection symptoms of seedlings are associated with the ability of the host to produce resin and the development of lignified zones (Cleary et al., 2012;Rishbeth et al., 1985), while the infection of excised shoots may represent the ability of the host to contain the mycelial propagation within the woody tissues (Heinzelmann & Rigling, 2016;Underwood & Pearce, 1992 (Adamson et al., 2015;Brodde et al., 2019). Because some of the largest and most productive maritime pine forests are located in the Atlantic distribution of the species (e.g., the Landes region in southwestern France), our results suggest that the expected increase of drought events in these populations will most likely cause severe damages due to their high susceptibility to D. sapinea. In the case of A. ostoyae, the main threat resides in the condition of the host. As mentioned before, a weaker host will be more susceptible to the fungus, and future extreme weather events are bound to weaken trees, also increasing pathogenicity of A. ostoyae (Kubiak et al., 2017). A mathematical model predicted a drastic northward shift of A. ostoyae in the northwestern United States for the years 2061-2080, leading to increased mortality of stressed and maladapted trees (Hanna et al., 2016). As suggested by this study, trees maladapted to new temperatures are also expected to be more susceptible to biotic stress. An interesting future avenue for further research will be the evaluation of susceptibility to both pathogens in the same season and ideally simultaneously on the same plant material to understand trade-offs in disease resistance.

| Genotype-phenotype association and candidate genes
By using a single-locus, two-step approach, we revealed significantly associated loci for all heritable traits under study. However, genotype effects were small, pointing to a highly polygenic nature of the studied traits, as often reported for adaptive traits in forest trees. In addition, for susceptibility to D. sapinea and A. ostoyae, no resistance alleles with major effects were detected. Nevertheless, information on genotype-phenotype associations can make genome-based phenotypic predictions more efficient by providing a set of SNPs that are more closely linked to causal SNPs than those selected just based on genome location (e.g., Westbrook et al., 2013). We retrieved annotations from Plomion et al. (2016) and found four non-synonymous SNPs significantly associated with duration of bud burst in 2015 (one locus), bud burst in 2017 (two loci) and needle discoloration caused by D. sapinea (one locus) under the additive genetic model (see Table 5).
The potential function of these genes has to be interpreted with caution as this information usually derives from studies in distantly related model species (typically Arabidopsis thaliana). Nevertheless, the locus (BX679001-1418), which was significantly associated with needle discoloration caused by D. sapinea, possibly codes for a translation initiation factor, eIF-5, that has previously been reported to be involved in pathogen-induced cell death and development of disease symptoms in A. thaliana (Hopkins et al., 2008). This gene deserves further attention in future studies addressing the genetic control of adaptive traits in conifers. Also, differential gene expression analyses could shed light on the genes and pathways involved in pathogen resistance in maritime pine, when applied to inoculated and control trees growing under controlled conditions in greenhouses or common gardens.
Based on a well-replicated clonal common garden and state-ofthe-art genotyping technology, we were able to study key adaptive traits in maritime pine and found evidence for non-synonymous mutations underlying genetic variation for some of these traits.
Association studies for highly polygenic traits are still challenging.

| CON CLUS IONS
In our study, we evaluated phenotypic variability and broad-sense heritability for height, spring phenology, pathogen susceptibility (as indicated by necrosis length and needle discoloration on excised branches), and incidence of processionary moth nests in a range-wide clonal common garden of maritime pine. We revealed strong genetic divergence at several adaptive traits, especially height and necrosis length caused by D. sapinea, across maritime pine populations. We have also shown that several adaptive traits in maritime pine were genetically correlated and that population variation was often established along climatic clines, in particular with maximum summer temperatures. The evolution of suits of functional traits along environmental clines is a common pattern in nature (e.g., Chapin et al., 1993;Reich et al., 1996). Currently, locally adapted populations are challenged by changing climate conditions, and emergent pests and pathogens expanding their range (Seidl et al., 2017). Susceptibility to D. sapinea was highest in the Atlantic maritime pine populations where it is expected to cause severe outbreaks due to increased incidence of drought events in future. Correlated selection for increased susceptibility needs to be considered in breeding programs aiming at increasing height growth in this species. In addition, opposing trends in pathogen susceptibility across maritime pine populations, for example, for D. sapinea and A. ostoyae (this study), and for the invasive pathogen F. circinatum (Elvira-Recuenco et al., 2014), challenge forest tree breeding and natural forest resilience. An improved understanding of integrated phenotypes, including responses to known pests and pathogens, and their underlying genetic architecture is fundamental to assist new-generation tree breeding and the conservation of valuable genotypes. Coupled with early detection methods (see, e.g., Kenis et al., 2018), knowledge on genetic responses to emerging pests and pathogens will help to ensure the health of forests in future. Finally, given the recent development of efficient technologies to combine functional genetic variants (e.g., genome editing), the identification of large numbers of polymorphisms associated with commercial traits is expected to contribute to a new era in plant breeding (Breeding 4.0, see Wallace et al., 2018).

ACK N OWLED G M ENTS
We

CO N FLI C T O F I NTE R E S T
We declare no conflict of interest.

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 are openly available in ["DRYAD"]