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

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.

INTRODUCTION

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., 2007, 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 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. Especially for model species, such as Populus sp. and Eucalyptus sp., association studies and complementary approaches revealed numerous genes involved in pathogen and pest resistance (e.g., Mhoswa et al., 2020; Muchero et al., 2018; Resende et al., 2017; Zhang et al., 2018).

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/ressources/memento‐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 extinction–recolonization events associated with past climatic changes as its host (Labbé et al., 2017). 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.

MATERIALS AND METHODS

Plant material and common garden measurements

A clonal common garden (CLONAPIN) was planted in 2011 in Bordeaux (44°44′43″N, 0°47′08″W), southwestern France (for details, see Rodríguez‐Quilón, 2017). It includes trees from 36 populations of maritime pine covering the whole natural species distribution (see Figure 1, Table S1.1, for number of individuals and genotypes, and population coordinates of 33 populations included in this study), representing all known differentiated gene pools in the species (Central Spain, southeastern Spain, Atlantic Iberian Peninsula, Atlantic France, Corsica, and Morocco; see 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. 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).

FIGURE 1

Map showing the origin of the Pinus pinaster populations included in the clonal common garden (part of CLONAPIN) established in Bordeaux, France (marked with a star). The coloring indicates the assignment of the populations to six neutral gene pools characterized by Jaramillo‐Correa et al. (2015)

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). Tree height was measured in 2015, four years after the establishment of the trial. On 11 weekly dates between April 3 and June 18, 2015, and on 13 weekly dates between March 23rd March and June 15th 2017, bud burst stage was evaluated using a phenological scale ranging from 0 to 5 (0: bud without elongation during winter, 1: elongation of the bud, 2: emergence of brachyblasts, 3: brachyblasts 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, mycelium‐side 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) for the 11 non‐Spanish and the 22 Spanish populations, respectively. Climate variables included monthly mean, highest and lowest temperatures, and mean monthly precipitation. Gonzalo’s (2007) model was favored for climate data in Spain because it considers a much denser network of meteorological stations than Worldclim, which is known to underperform, particularly for precipitation estimates, in this region (see Jaramillo‐Correa et al., 2015).

DNA extraction and SNP genotyping

Needles were collected from one replicate per genotype (N = 416, including all genotypes used for pathogen inoculations) and desiccated using silica gel. Genomic DNA was extracted using the Invisorb® DNA Plant HTS 96 Kit/C kit (Invitek GmbH, Berlin, Germany). An Illumina Infinium SNP array developed by Plomion et al. (2016) was used for genotyping. Apart from potentially neutral genetic polymorphisms, this array comprises SNPs from candidate genes that showed signatures of natural selection (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 (Plomion et al., 2016). 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, µ denotes the overall phenotypic mean, block represents the fixed effect of experimental block i, gp is the random effect of gene pool g, gp(pop) is the random effect of pop j nested within gene pool g, gp(pop(genotype)) 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 using Markov chain Monte Carlo (MCMC) methods implemented in the R package MCMCglmm (Hadfield, 2010) using R v.3.4.1 (R Development Core Team, 2017). All analyzed traits presented a Gaussian distribution with the exception of presence of pine 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 (108) 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 is the variance among genotypes within populations and gene pools and the residual variance. In addition, we computed the variance ratios associated with the population () and gene pool () effects as: where is the variance among populations within gene pools, and 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 is the overall variance among populations, and 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 correlations between traits that can constrain adaptive evolution, Etterson & Shaw, 2001). Pearson's correlation coefficient was also used to characterize correlations of trait variation across populations, that is, considering the combined genotype and population effects. These correlations reflect both the effect of pleiotropy and correlation of causative loci within‐populations, on one hand, and the effect of distinct evolutionary histories experienced by populations over the range of the species distribution, on the other hand. BLUPs were obtained from the model including genotype and population effects but not the gene pool effects (second set, equations 3 and 4), as they produced more precise estimates. In addition, we employed linear mixed models to test the effect of environmental variables on the population BLUPs while taking into account the assignment of each population to one of the six gene pools as random factor. First, we tested each model against a null model without fixed effects to identify the significant associations. Then, for significant environmental associations, we calculated R 2 as the proportion of variance explained by the fixed predictors using Nakagawa & Schielzeth's approach (Johnson, 2014; Nakagawa & Schielzeth, 2013) implemented in the function “r2beta” of the r2glmm package (Jaeger, 2017).

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 p‐value 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 approach (Bayesian Association with Missing Data, BAMD; Li et al., 2012; Quesada et al., 2010) in R to estimate single‐locus genotype effects under three genetic models accounting for additive, overdominance and dominance effects, respectively (as in 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.

RESULTS

Phenotypic variability, broad‐sense heritability, and genetic differentiation

Broad‐sense heritability was highest for bud burst in 2015 (H 2: 0.310, CI [0.272–0.358]), intermediate for height (H 2: 0.219 [0.186–0.261]) and lowest for necrosis length caused by A. ostoyae (H 2: 0.021 [0.004–0.121], Table 1). Broad‐sense heritability of susceptibility to D. sapinea, assessed as the necrosis length, was also significant (H 2: 0.152 [0.038–0.292]), with trees from northern Africa and southern Spain showing shorter necrosis length than trees from Atlantic populations (Figure 2a). Broad‐sense heritability of needle discoloration caused by D. sapinea had similar magnitude, but was not significant (H 2: 0.123 [0.000–0.233]). Necrosis length caused by A. ostoyae was most pronounced in southern populations, especially in Morocco and southern Spain, while less pathogen growth was observed in northern populations, for example, in St. Jean de Monts from the French Atlantic gene pool (Figure 2b). Incidence of pine processionary moth nests in the common garden was not heritable (H 2: 0.001 [0.000–0.207]). The tallest trees were found in populations from Atlantic France and Atlantic Iberian Peninsula gene pools and for one of the Corsican populations (Pinia), whereas the shortest ones originated from southeastern Spain and Morocco (Figure S5.1).

TABLE 1

Broad‐sense heritability, variance ratios, and genetic differentiation of adaptive traits in Pinus pinaster

	Phenotypic mean	Variability	H 2	H 2 p	H 2 gp	Q ST
height (cm)	170.647	± 48.228	0.219 [0.186–0.261]	0.069 [0.028–0.132]	0.254 [0.113–0.714]	0.549 [0.392–0.662]
bb2015 (dd)	1311.95	± 82.411	0.310 [0.272–0.358]	0.058 [0.022–0.121]	0.194 [0.06–0.581]	0.275 [0.186–0.443]
dbb2015 (dd)	814.713	± 116.708	0.165 [0.099–0.212]	0.040 [0.000–0.111]	0.000 [0.000–0.287]	0.191 [0.100–0.332]
bb2017 (dd)	1286.245	± 79.853	0.230 [0.181–0.261]	0.031 [0.010–0.088]	0.066 [0.000–0.302]	0.191 [0.106–0.404]
dbb2017 (dd)	901.149	± 78.922	0.237 [0.195–0.297]	0.000 [0.000–0.063]	0.231 [0.102–0.676]	0.463 [0.293–0.579]
A. ostoyae necrosis (mm)	48.533	± 29.625	0.021 [0.004, 0.121]	0.066 [0.018, 0.203]	0.066 [0.018, 0.203]	0.217 [0.041–0.787]
D. sapinea necrosis (mm)	43.348	± 17.931	0.152 [0.038–0.292]	0.000 [0.000–0.282]	0.002 [0.000–0.602]	0.636 [0.349–1.000]
D. sapinea disc	0 ‐ no disc: 183 1 ‐ low: 123 2 ‐ medium: 141 3 ‐ high: 9		0.123 [0.000–0.233]	0.001 [0.000–0.145]	0.002 [0.000–0.296]	0.093 [0.000–0.752]
Processionary	1 ‐ presence: 48 0 ‐ absence: 3282		0.001 [0.000–0.207]	0.001 [0.000–0.128]	−0.004 [0.000–0.380]	0.006 [0.000–0.985]

Variability refers to the standard deviation of the raw phenotypic data. H 2, broad‐sense heritability; H 2 p, variance ratio of the population effect (see main text for ratio definition); H 2 gp, variance ratio of the gene pool effect (see main text for ratio definition); Q ST, population differentiation at quantitative traits.

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.

FIGURE 2

Stripchart of combined genotype and population BLUPs (best linear unbiased predictors) for necrosis length caused by D. sapinea (a) and A. ostoyae (b). Populations were assigned to one of the six gene pools identified by Jaramillo‐Correa et al. (2015), corresponding to the six colors in the figure, and ordered by latitude (north to south). Symbols indicate different populations of the same gene pool. Black lines indicate the average necrosis length in each population

Most traits showed strong differences among gene pools, although the variance ratio associated with the gene pool effect was characterized by wide confidence intervals (Table 1). Based on overall population genetic differentiation, we observed significant Q ST values (ranging from 0.191 for bud burst in 2017 and duration of bud burst in 2015, to 0.636 for necrosis length caused by D. sapinea) indicating strong population differentiation (Table 1). Global F ST calculated using the available SNPs was 0.109 ([0.013; 0.325], p‐value < 0.001) which is significantly lower than the Q ST estimates obtained for height and necrosis length caused by D. sapinea (Table 1).

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 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, p‐value < 0.001), and positive for necrosis length caused by D. sapinea (0.679, p‐value < 0.001).

TABLE 2

Pearson's correlation coefficients of the best linear unbiased predictors (BLUPs) of genotype (i.e., genetic correlations, above diagonal) and combined population and genotype effects (below diagonal)

	Height	bb2015	dbb2015	bb2017	dbb2017	D. sapinea necrosis	D. sapinea disc	A. ostoyae necrosis
Height		−0.172**	−0.132*	−0.178**	0.439***	0.111	0.082	−0.010
bb2015	0.194***		0.798***	0.392***	−0.292***	0.118	0.041	−0.017
dbb2015	0.117*	0.846***		0.380***	−0.233***	0.061	0.065	−0.110
bb2017	0.021	0.593***	0.547***		−0.154**	−0.023	0.069	−0.092
dbb2017	0.800***	0.147**	0.079	0.018		−0.032	−0.093	−0.030
D. sapinea. necrosis	0.679***	0.564***	0.511***	0.309***	0.639***		0.231*	−0.009
D. sapinea disc	0.414***	0.470***	0.492***	0.337***	0.392***	0.727***		−0.149
A. ostoyae necrosis	−0.653***	−0.321*	−0.284*	−0.056	−0.530***	−0.692***	−0.367**

Abbreviations: bb, bud burst; dbb, duration of bud burst; disc, needle discoloration.

Significance levels after false discovery rate (FDR) correction: *< 0.05; **< 0.01; ***0.001.

FIGURE 3

Correlation of trait variation across Pinus pinaster populations based on BLUPs (best linear unbiased predictors) for necrosis length caused by Diplodia sapinea and Armillaria ostoyae. A linear trend line is also shown (Pearson's correlation coefficient = −0.692, 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). The strongest relationship was found for height and maximum temperatures in July indicating a negative correlation (R 2 = 0.709 [0.556–0.831], p‐value <0.001) while the variance of, for example, needle discoloration caused by D. sapinea explained by maximum temperatures in July was lower but also negative (R 2 = 0.524 [0.128–0.841], p‐value = 0.004; Table 3, Figure 4). Also, geographical variables, such as altitude, latitude, and longitude, showed significant correlation with all phenotypes. Interestingly, necrosis length caused by A. ostoyae showed a latitudinal cline (R 2 = 0.750, p‐value = 0.016) but only a marginally non‐significant positive correlation with maximum temperatures in July (R 2 = 0.514, p‐value = 0.066).

TABLE 3

Correlations between adaptive phenotypes and environmental variables, as identified using linear mixed‐effects models that account for gene pool as a random factor

Trait	Environmental factor	Direction of correlation	R 2	p
Height	Altitude	−	0.650	<0.001
	Tmax July	−	0.709	<0.001
	Tmax August	−	0.201	0.023
bb2015	Altitude	−	0.200	0.020
dbb2015	Longitude	−	0.245	0.044
bb2017	Longitude	−	0.262	0.028
dbb2017	Latitude	+	0.290	0.042
	Altitude	−	0.430	0.001
	Tmax July	−	0.212	0.044
	Tmax August	−	0.193	0.033
	Prec in August	+	0.123	0.027
	Prec in September	+	0.123	0.027
D. sapinea necrosis	Latitude	+	0.341	0.031
	Altitude	−	0.482	0.018
	Tmax July	−	0.465	0.034
	Tmax August	−	0.394	0.036
D. sapinea disc	latitude	+	0.284	0.013
	Tmax June	−	0.310	0.009
	Tmax July	−	0.524	0.004
	Tmax August	−	0.462	0.003
	Prec. in September	−	0.315	0.004
A. ostoyae necrosis	Latitude	−	0.750	0.016
	Altitude	−	0.614	0.040

Only models that performed significantly better than a null model excluding the fixed effect are reported. R 2, the proportion of variance explained by the fixed predictor in the mixed‐effects model. For necrosis length caused by A. ostoyae, a simple linear regression model was used as each gene pool was only represented by a single population.

Abbreviations: bb, bud burst; dbb, duration of bud burst; disc, needle discoloration; Tmax, maximum temperature; prec, precipitation.

FIGURE 4

Stripchart of combined genotype and population BLUPs (best linear unbiased predictors) for necrosis length caused by D. sapinea plotted against the maximum temperature in July in the population of origin. Populations were assigned to one of the six gene pools identified by Jaramillo‐Correa et al. (2015), corresponding to the six colors in the figure. Symbols indicate different populations of the same gene pool. Black lines indicate the average necrosis length in each population

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 non‐coding 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).

TABLE 4

Single nucleotide polymorphisms (SNPs) significantly associated with height, spring phenology, and pathogen susceptibility traits under the additive genetic model, as identified by a two‐step approach based on mixed‐effects linear models (MLMs) implemented in Tassel and the Bayesian framework in BAMD (BMLMs)

Trait	SNP name	SNP motif	Site annotation	LG	MAF	MLMs			BMLMs
						F	p	R 2	Mean	Lower 95% CI	Upper 95% CI
Height	BX249583‐420	[A/G]	unk	1	0.266	7.383	0.001	0.027	3.956	1.162	6.762
	BX251999‐509	[A/T]	unk	9	0.296	6.576	0.002	0.024	2.749	0.262	5.255
	CT2714‐442	[T/A]	unk		0.381	5.683	0.004	0.021	−2.920	−5.335	−0.503
	F51TW9001AZG2W−933	[C/G]	unk	4	0.438	9.991	<0.001	0.036	3.480	0.931	6.083
	sp_v3.0_unigene17345‐1191	[T/G]	nc	9	0.344	5.977	0.003	0.022	3.232	0.588	5.871
bb2015	BX249218‐322	[A/C]	nc		0.315	6.499	0.002	0.03	7.230	2.107	12.440
	BX249671_307	[T/C]	unk	7	0.397	6.438	0.002	0.029	6.129	0.763	11.450
	BX253890‐151	[T/G]	nc	12	0.157	6.350	0.002	0.029	10.95	4.298	17.552
	CL2033CT1302CN1398‐513	[A/G]	nc	1	0.408	5.436	0.005	0.025	7.134	2.001	12.321
	CL544Contig1_03. Pipn−84	[T/G]	unk		0.135	9.107	<0.001	0.041	10.18	2.922	17.579
	FN692276‐550	[A/G]	unk	12	0.402	7.903	<0.001	0.036	5.667	0.664	10.774
	i13066s710	[A/C]	nc		0.242	6.107	0.002	0.028	8.365	2.328	14.407
	i16267s380	[A/G]	unk	2	0.411	7.388	0.001	0.033	−10.55	−15.885	−5.168
	LP3‐3–298	[C/G]	unk		0.143	5.538	0.004	0.025	8.510	1.696	15.302
bb2017	0_12730_01_contig1‐159	[A/C]	unk	12	0.379	9.749	<0.001	0.044	4.003	1.509	6.499
	AL750545‐695	[T/A]	non‐syn	1	0.487	5.754	0.003	0.026	4.131	1.375	6.800
	AL750773_910	[T/A]	unk	3	0.499	5.580	0.004	0.026	−4.048	−6.810	−1.339
	BX249816‐2143	[A/G]	non‐syn	7	0.269	5.783	0.003	0.026	5.787	2.792	8.839
	CT576106‐142	[C/G]	unk	1	0.180	7.147	0.001	0.032	5.590	2.165	9.029
	CT577489‐1569	[A/C]	unk		0.203	5.826	0.003	0.026	−4.606	−7.559	−1.602
	F7JJN6E01B7BCW−157	[A/G]	syn	5	0.117	5.591	0.004	0.025	7.708	3.801	11.533
	FM945796‐840	[T/G]	unk		0.213	6.284	0.002	0.028	−4.613	−7.748	−1.473
	i10996s1211	[T/C]	unk		0.301	8.281	<0.001	0.037	4.295	1.467	7.149
dbb2015	CL2507CT3369CN3610‐193	[T/A]	nc	1	0.387	5.406	0.005	0.025	8.981	2.203	15.730
	F51TW9001BWV4H−219	[T/C]	non‐syn		0.462	5.941	0.003	0.028	7.101	0.458	13.72
	F51TW9002FPGRE−170	[A/G]	nc		0.346	5.510	0.004	0.026	9.467	2.507	16.580
dbb2017	AL749850_679	[A/G]	unk		0.402	5.449	0.005	0.021	−4.992	−9.318	−0.654
	BX681656‐960	[A/C]	syn	5	0.481	8.050	<0.001	0.030	−6.416	−10.154	−2.669
	CT580064‐331	[A/G]	unk		0.442	6.236	0.002	0.024	−4.363	−7.934	−0.822
	CT582680‐451	[A/C]	unk		0.201	7.383	0.001	0.028	−10.48	−15.541	−5.470
	F51TW9001BAW7V−405	[A/G]	unk	12	0.163	9.088	<0.001	0.034	7.287	2.665	11.854
	i17647s350pg	[G/C]	unk		0.157	5.847	0.003	0.022	5.747	0.878	10.611
D. sapinea disc	BX679001‐1418	[T/C]	non‐syn	7	0.192	5.544	0.005	0.049	−0.054	−0.102	−0.006
	CT574726_30	[T/C]	unk		0.435	5.539	0.005	0.049	−0.051	−0.089	−0.014
	i11062s233	[G/C]	syn	7	0.293	5.596	0.005	0.050	−0.048	−0.090	−0.007
D. sapinea necrosis length	AL750492_595	[T/C]	unk	1	0.121	5.887	0.004	0.033	0.897	0.081	1.713
	F51TW9002FT2ZF−1060	[T/G]	unk	12	0.485	8.614	<0.001	0.048	−1.279	−1.992	−0.574
A. ostoyae necrosis length	F51TW9001AI9YZ_1847	[A/G]	unk	7	0.273	5.928	0.005	0.081	−0.735	−1.345	−0.129
	F51TW9001CXU1D_1264	[T/C]	unk	6	0.364	6.594	0.003	0.09	−0.997	−1.787	−0.225

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.

TABLE 5

Annotation for SNPs significantly associated with adaptive traits under the additive model and coding for a non‐synonymous amino acid change, as retrieved from Plomion et al. (2016)

Trait	SNP name	Motif	Protein change	Putative protein function
dbb2015	F51TW9001BWV4H−219	[T/C]	Asparagine ‐ Serine	LANC‐like domain containing protein
bb2017	BX249816‐2143	[A/G]	Glutamate/Glutamine ‐ Arginine	putative Cdc2‐related protein kinase CRK2
bb2017	AL750545‐695	[A/T]	Glutamate/Glutamine ‐ Valine	Catalase
D. sapinea disc.	BX679001‐1418	[T/C]	Isoleucine ‐ Valine	Translation initiation factor eIF−5

Abbreviations: bb, bud burst; dbb, duration of bud burst; disc, needle discoloration.

FIGURE 5

Density plots (left) of the effect sizes based on 20,000 BAMD simulations with 95% confidence intervals (dashed lines) and genotypic effects (box plots, right) for four non‐synonymous single nucleotide polymorphisms (SNPs) showing significant association with (from top to bottom) needle discoloration caused by Diplodia sapinea, duration of bud burst in 2015, and bud burst in 2017 (two SNPs) in Pinus pinaster

DISCUSSION

In the current context of climate change, understanding the genetic basis of adaptive traits in tree species is key for an informed forest management and the development of new breeding applications. In this study, we assessed variation in maritime pine for height and spring phenology traits, as well as for incidence of pine processionary moth nests and response to two pathogenic fungi, D. sapinea and A. ostoyae, at the range‐wide scale. Novel inoculation protocols were developed for D. sapinea and A. ostoyae based on excised branches and applied to trees grown in a clonal common garden. For the first time, we estimated broad‐sense heritability of pine susceptibility for these two major fungal pathogens (evaluated as necrosis length on detached branches), as well as the variance ratios associated with the population and gene pool effects. Most of these adaptive traits showed significant broad‐sense heritability. Variation of necrosis length caused by pathogens among geographical provenances, as well as duration of bud burst in 2017, followed a latitudinal gradient, possibly corresponding to a climatic gradient, but in opposite directions for the two pathogens. A genetic association approach revealed several loci significantly associated with height, spring phenology, and traits related to pathogen infection in maritime pine. This information will be useful for current efforts to implement breeding strategies based on genome prediction (i.e., genomic selection) in this species. The presence of pine processionary moth nests evaluated in the common garden was not heritable but future studies should consider the level of infestation and/or assess damage more precisely.

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 the population level, for example, drift or different levels of gene flow. Nevertheless, they provide important insights on trait variation at the range‐wide geographical scale. The variance ratio associated with the gene pool effect had very wide confidence intervals indicating that highly differentiated units were mixed at this level. Already, Rodríguez‐Quilón et al. (2016) described strong quantitative genetic differentiation for height and survival within gene pools of maritime pine. While gene pools may reflect 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 moderate broad‐sense heritability of 0.219 for height, well in line with estimates in other conifer species, for example, ranging from 0.21 in Pinus taeda to 0.78 in Picea abies (reviewed in Lind et al., 2018) and from 0.148 to 0.282 in maritime pine saplings depending on the common garden site and the provenance (Rodríguez‐Quilón et al., 2016). We found significant genetic correlation between height and spring phenology indicating possibly pleiotropic effects of genes involved in the expression of these traits. The direction of genetic correlations and of trait variation across populations showed opposing trends in some cases. Height showed negative genetic correlations with spring phenology (except for duration of bud burst 2017) but positive correlations with the same traits when considering 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 (Alía et al., 2014; Jaramillo‐Correa et al., 2015), and has thus been used in combination with genetic markers to identify relevant conservation units in maritime pine (Rodríguez‐Quilón et al., 2016). 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 (2015 or 2017). Differentiation (Q ST) for bud burst varied from 0.191 to 0.275, which is comparable to a mean of 0.249 for bud flush averaged over several forest tree species (reviewed in Alberto et al., 2013). In our study, trees originating from northern populations flushed later than trees from southern populations. Similar clines have been observed, for example, in Castanea sativa, another south European tree species (Míguez‐Soto et al., 2019), and in French provenances of Quercus petraea (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.

Broad‐sense heritability for susceptibility to D. sapinea based on necrosis length (H 2 = 0.152) was much lower than, for example, narrow‐sense heritability of lesion length caused by Fusarium circinatum (h 2 = 0.618, Elvira‐Recuenco et al., 2014) in P. pinaster. However, this low but significant heritability indicated the potential for an evolutionary response of maritime pine to D. sapinea and that inoculations on excised branches were suitable to capture a genetic component of this trait. Moreover, Q ST was highest for necrosis length caused by this pathogen (QST = 0.636) in our study, which was also significantly higher than neutral‐marker F ST. This result pointed to adaptive divergence and substantial variability in susceptibility to D. sapinea across the geographical range of maritime pine. This population variability could be harnessed in maritime pine breeding programs to increase resistance to the disease (Bouffier et al., 2008; Ingvarsson & Dahlberg, 2019). Remarkably, we observed a negative correlation of trait variation across populations for necrosis length caused by D. sapinea and A. ostoyae. This direct comparison should be taken with caution as the inoculation experiments were carried out at different time points (D. sapinea in summer and A. ostoyae in autumn), which is likely to affect the host response (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). However, correlation with climate across populations for susceptibility to both D. sapinea and A. ostoyae (as well as for other adaptive traits such as height and spring phenology) indicated similar environmental clines driving differentiation at these traits. Especially, maximum temperatures during summer months showed significant correlations with genetic variability of susceptibility to pathogens across maritime pine populations. Trees from populations with high maximum summer temperatures were less susceptible to D. sapinea. This result can be interpreted in different ways: (1) If we assume that D. sapinea is native in Europe and northern Africa, the pathogen pressure can be expected to be stronger in southern regions, with a climate more favorable to D. sapinea pathogenic outbreaks, triggered by stress in the host plant, especially by droughts (Luchi et al., 2014). Maritime pine populations growing in these regions (e.g., Morocco and southern Spain) would then be more likely to have evolved resistance to the disease. In contrast, trees from populations where severe drought periods have most likely not been common so far (e.g., Atlantic populations from the Iberian Peninsula and France) would be more susceptible. (2) If maritime pine and D. sapinea have not had sufficient time to co‐evolve (e.g., if D. sapinea was not native to Europe and northern Africa as hypothesized by Adamson 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; de la Mata et al., 2014; Gaspar et al., 2013). 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) showed that the metabolomes in needles of maritime pine trees from populations with distinct geographic origin (notably Atlantic versus Mediterranean provenances) were quite differentiated, with flavonoids showing a significant correlation with the water regime of the population of origin. However, the expression of metabolites is organ specific (de Miguel et al., 2016) and knowledge about secondary metabolites involved in resistance to D. sapinea is still lacking.

A study on the invasive pathogen Fusarium circinatum, which did certainly not co‐evolve with maritime pine, also revealed a geographic cline in susceptibility, with Atlantic maritime pine populations showing less susceptibility than Moroccan populations (Elvira‐Recuenco et al., 2014). A similar pattern was observed for A. ostoyae in our study. Our results indicated that maritime pine from southwestern France, where A. ostoyae outbreaks have been frequently reported (Labbé et al., 2015), may have developed some resistance or might show exapted resistance to the disease. Considering the absence of reports of A. ostoyae from the south of the Iberian Peninsula (Marxmüller & Guillaumin, 2005), which is in line with the species’ preference for humid forest sites (Cruickshank et al., 1997; Heinzelmann et al., 2019), trees in Morocco and southern Spain have most likely never co‐evolved with this pathogen. However, a study by Guillaumin et al. (2005) on the mortality of potted maritime pine plants revealed an opposite pattern to the one found in our study, with the Landes population in Atlantic France being the most susceptible and the Moroccan population the least susceptible to A. ostoyae. Moreover, Zas et al. (2007) found much higher (narrow‐sense) heritability (h 2 = 0.35) for mortality due to A. ostoyae in an infected progeny trial of maritime pine seedlings than the values of broad‐sense heritability of necrosis length found in our study. These contrasting results call for further research to improve our understanding of the different components that affect resistance against A. ostoyae and which might vary among host populations. Armillaria ostoyae is a root pathogen and a critical point during natural infection that could be key for resistance mechanisms is the penetration into the root (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). In future, it would therefore be important to carry out root inoculations on potted seedlings or young trees from range‐wide maritime pine populations to confirm the patterns observed in this study. Another relevant factor to consider is the particular common garden environment, as it provided similar environmental conditions compared to their place of origin for some of the populations, especially the populations from southwestern France which might be locally adapted, but very different to those of other populations (e.g., those originating in Morocco). This could have affected trait variation due to genotype by environment (G x E) interactions. The expression of trait genetic differences dependent upon trial sites can, for example, lead to reduced heritability or extreme phenotypes only visible under certain environmental conditions or in certain years. Comparing results from several trial locations with contrasting environmental conditions covering the natural variability of the species would therefore be important in future studies.

Suitable strategies to evaluate susceptibility to D. sapinea and A. ostoyae will become increasingly important as climate change increases pathogen pressure. Droughts are expected to become more frequent throughout Europe (IPCC, 2014), which will most likely trigger D. sapinea outbreaks also in regions where the pathogen has not caused severe disease symptoms so far. Recently, a northward expansion of D. sapinea outbreaks in Europe, probably driven by higher spring temperatures, has been reported; these outbreaks have caused severe damage on P. sylvestris in Sweden and eastern Baltic countries (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‐of‐the‐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. Lind et al. (2017) reported an average of 236 SNPs (out of 116,231 tested, i.e., 0.203%) associated with each of four fitness‐related traits in Pinus albicaulis by detecting signals of significantly higher covariance of allele frequencies. In the near future, multilocus association methods should be used to reveal genome‐wide loci with non‐zero effects for polygenic traits in forest trees (De la Torre et al., 2019; Goldfarb et al., 2013).

CONCLUSIONS

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).

CONFLICT OF INTEREST

We declare no conflict of interest.

Supplementary Material

Click here for additional data file.

Table S6

Click here for additional data file.