Ecological and evolutionary drivers of phenotypic and genetic variation in the European crabapple [Malus sylvestris (L.) Mill.], a wild relative of the cultivated apple

Abstract

Background and Aims

Studying the relationship between phenotypic and genetic variation in populations distributed across environmental gradients can help us to understand the ecological and evolutionary processes involved in population divergence. We investigated the patterns of genetic and phenotypic diversity in the European crabapple, Malus sylvestris, a wild relative of the cultivated apple (Malus domestica) that occurs naturally across Europe in areas subjected to different climatic conditions, to test for divergence among populations.

Methods

Growth rates and traits related to carbon uptake in seedlings collected across Europe were measured in controlled conditions and associated with the genetic status of the seedlings, which was assessed using 13 microsatellite loci and the Bayesian clustering method. Isolation-by-distance, isolation-by-climate and isolation-by-adaptation patterns, which can explain genetic and phenotypic differentiation among M. sylvestris populations, were also tested.

Key Results

A total of 11.6 % of seedlings were introgressed by M. domestica, indicating that crop–wild gene flow is ongoing in Europe. The remaining seedlings (88.4 %) belonged to seven M. sylvestris populations. Significant phenotypic trait variation among M. sylvestris populations was observed. We did not observe significant isolation by adaptation; however, the significant association between genetic variation and the climate during the Last Glacial Maximum suggests that there has been local adaptation of M. sylvestris to past climates.

Conclusions

This study provides insight into the phenotypic and genetic differentiation among populations of a wild relative of the cultivated apple. This might help us to make better use of its diversity and provide options for mitigating the impact of climate change on the cultivated apple through breeding.

Population structure, isolation by distance, isolation by ecology, local adaptation, climate change, apple tree, crop wild relatives

INTRODUCTION

Knowledge of the spatial phenotypic and genetic variation among populations is essential for understanding the ecological (biotic and abiotic factors) and evolutionary (gene flow, selection, drift and mutation) processes involved in population divergence and adaptation (Savolainen et al., 2013; Sork, 2018). Plant species distributed across climatic gradients typically experience spatial variation in selection, genetic drift and gene flow, which drives genetic and phenotypic divergences among populations (Svenning et al., 2015). Climate influences demographic processes such as population expansion and contraction, the extent of gene flow among populations and, ultimately, the degree of genetic divergence among populations (Edwards et al., 2022). For instance, changes in the climate since the Last Glacial Maximum (LGM) 20 000 years ago have driven the genetic composition of the European crabapple and many other tree species (Comes and Kadereit, 1998; Kremer et al., 2002; Pyhäjärvi et al., 2008; Cornille et al., 2013a; Riordan et al., 2016; Lander et al., 2021; Parisod, 2021; Yamada et al., 2021). Climate can also shape phenotypic variation among populations. Populations occurring in the same climate can share physiological tolerances to climatic conditions, including plant carbon uptake via photosynthesis. Carbon uptake traits condition plant size and growth, reproduction and survival in different climatic conditions (Hutyra et al., 2007; Nicotra et al., 2010; Pasho et al., 2012; Way and Montgomery, 2015; Mercado et al., 2018; Hartmann et al., 2020; Kühn et al., 2021). Variation in plant carbon uptake in response to climatic conditions can result from phenotypic plasticity, i.e. the ability of individual genotypes to produce different phenotypes when exposed to other environmental conditions, in this case, climate (Dusenge et al., 2019). In some cases, local climate can impose divergent selection on carbon uptake traits, which leads to reproductive isolation among populations; loci under adaptive divergence act as a local barrier to gene flow (Keller et al., 2011; Franks et al., 2014; Aitken and Bemmels, 2015; Ramírez-Valiente et al., 2017; Alexandre et al., 2020). Climate can therefore lead to a long-term reduction in gene flow and local adaptation. Whether the phenotypic variation observed in species distributed across large climatic ranges results from their demographic and/or adaptive histories remains an intense topic of investigation (Li et al., 2012; Tiffin and Ross-Ibarra, 2014; Briscoe Runquist et al., 2020). Investigating this question can help us to predict how plants might respond to climate change and how species adapt to their environment.

There are multiple ways to investigate whether genetic and phenotypic variation among populations distributed across climatic gradients results from selection, genetic drift and/or gene flow. A first step is to use a common garden experiment to investigate the genetic basis of phenotypic variation among populations. Different populations occurring across a climatic gradient might display clinal variation, i.e. differences in a trait that might result from phenotypic plasticity or local genetic adaptation (Savolainen et al., 2013; de Villemereuil et al., 2016). Measuring candidate traits for adaptation to climate [e.g. phenology (Brachi et al., 2013) or traits related to plant carbon uptake (Savolainen et al., 2013; de Villemereuil et al., 2016)] in individuals from different populations in the same environmental conditions can help us to elucidate the genetic basis of phenotypic variation across populations without the confounding effects of the environment. Ideally, common garden trials should include the main genetic groups across the species distribution (de Villemereuil et al., 2020). The association of neutral genotypic variation with phenotypic variation can also be used as evidence of adaptive divergence among populations (Shafer and Wolf, 2013; Wang and Bradburd, 2014). The correlation between neutral genetic differentiation and environmental or phenotypic divergence among populations, independent of geographical distance, referred to as isolation by ecology (IBE), is an extension of the isolation-by-distance (IBD) model (Wright, 1943). Patterns of IBE have been increasingly used as an indicator of adaptive divergence between populations (Shafer and Wolf, 2013). In the IBE model, natural selection, which results from several factors, including climate, can indirectly increase neutral genetic and phenotypic differentiation between populations by promoting general barriers to gene flow (Nosil et al., 2009; Orsini et al., 2013; Shafer and Wolf, 2013; Wang and Bradburd, 2014). Although it can be challenging to map the processes underlying IBE patterns, testing for it is valuable to understand better how natural selection shapes neutral genetic and phenotypic variation. Therefore, evidence from common garden experiments and IBE patterns can contribute to understanding how genotypes, phenotypes and the environment interact to influence population divergence and, potentially, local adaptation.

Fruit trees are a significant component of terrestrial ecosystems (Petit and Hampe, 2011). They are grown in managed plantations and orchards to provide a variety of economically important products (Boyd et al., 2013). Recent breeding efforts have involved the repeated use of a limited number of commercial cultivars, leading to a reduction in genetic diversity that can lead to the loss of valuable alleles of genes that are not directly targeted by human selection (Myles et al., 2011; Warschefsky and von Wettberg, 2019; Migicovsky et al., 2021). Wild relatives of crop fruit trees [crop wild relatives (CWRs)] harbour phenotypic and genetic diversity that is potentially of high valuable for future breeding programmes in the context of climate change (Zhang et al., 2017; Hoban et al., 2018; Hübner and Kantar, 2021). However, phenotypic variation in CWRs in relationship to climate variation has rarely been studied in fruit trees (Kremer and Hipp, 2019). Key traits to study in this context are related to plant carbon uptake. Climate impacts plant carbon uptake (Aubin et al., 2016), which impacts fruit quality characteristics and production (Demestihas et al., 2017). These are topical issues because native CWRs can be threatened by crop-to-wild gene flow from nearby domesticated trees (Delplancke et al., 2011; Cornille et al., 2015; Diez et al., 2015; Feurtey et al., 2017; Flowers et al., 2019; Liu et al., 2019). Therefore, studying the genetic and phenotypic variation among CWR fruit tree populations is timely to guide future breeding programmes. It might also contribute to our understanding of the evolutionary and ecological drivers of population divergence, such as climate.

The European crabapple, Malus sylvestris (L.) Mill., is a CWR of the cultivated apple, Malus domestica (Cornille et al., 2012, 2014, 2019; Peace et al., 2019). Substantial crop-to-wild gene flow has been observed across M. sylvestris populations in Europe [≤23.1 % of naturally occurring individuals are introgressed by M. domestica (Cornille et al., 2015)]. Crop–wild hybrids sampled in a forest in France and grown in controlled conditions showed higher growth rates than wild seedlings (Feurtey et al., 2017). Population genetic analyses also identified five pure (i.e. not introgressed by M. domestica) populations in Scandinavia, western France, eastern France, Eastern Europe and Italy (Cornille et al., 2015). These five populations resulted from past contractions and expansions associated with the LGM (Cornille et al., 2013a; Cornille et al., 2015). It is not known whether these five populations distributed across a large area with different climatic conditions show a phenotypic variation that could result from local adaptation to past and/or present climates.

We investigated the spatial phenotypic and genetic variation among populations of a wild contributor to the cultivated apple genome, M. sylvestris, to test for adaptive divergence, particularly to climate. Plant growth and traits related to carbon uptake were measured in 584 M. sylvestris seedlings grown in controlled conditions and genotyped using 13 microsatellite markers. We first assessed the genetic status of each seedling (pure vs. crop–wild hybrid). Then, we compared growth traits and traits related to plant carbon uptake among seedlings from different European genetic groups. We also formally tested the impact of geography (IBD) and ecology (IBE tested with phenotypic traits and climate) on genetic variation observed from 13 microsatellite markers. We investigated the following questions. First, do growth rates and carbon uptake traits vary between populations of the European crabapple, and are those traits heritable and therefore good candidates for the responses to selection? Second, is there any association between phenotypic variation and genetic variation, taking into account geographical distance, which would suggest adaptive divergence? Third, do we detect a pattern of isolation by climate in the European crabapple that could suggest local adaptive divergence to climate?

MATERIALS AND METHODS

Plant material, experimental design and trait measurements

A total of 584 seeds were collected from 90 M. sylvestris mother trees (between 3 and 15 seeds per mother tree; Supplementary data Table S1) from 22 geographical sites in Europe representing the five main genetic groups previously detected by Cornille et al. (2015): Austria (n = 89, two sites), Denmark (n = 91, three sites), Spain (n = 39, one site), France (n = 220, eight sites), Italy (n = 32, one site) and Romania (n = 117, seven sites) (Supplementary data Table S1). The length, width and weight of 30 seeds per mother tree were measured with an Opto-Agri (Optomachine, Riom, France).

In mid-April 2019, 584 seeds were washed, sterilized (in 0.5 % chlorine for 20 min), and stratified for 3 months at 4 °C in the dark in a mix of damp sand and vermiculite. Then, seeds were sown in jiffy pellets, and each pellet was placed randomly in a 20-hole array. Seeds were grown in controlled conditions for 2 months (from mid-July to mid-September 2019: 22 ± 1 °C, 60 ± 5 % relative humidity, a 16 h:8 h (light:dark) photoperiod, and a light level of 40–60 µmol m⁻² s⁻¹). Each 20-hole array was rotated daily in the growth chamber to avoid any micro-environmental variation in plant response, and plants were watered weekly.

During the 2-month experiment, the number of leaves and height of each seedling was recorded. Owing to the low germination rate, some accessions could not be evaluated (i.e. plant height and number of leaves could not be recorded for 19 of 584 seedlings, resulting in n = 565; Table 1). Seedlings were measured every 2–3 days, starting from days 7–11 after commencing the experiment.

Table 1.

Number of Malus sylvestris seedlings used in this study for population genetic analyses inferred with STRUCTURE for K = 8, with 13 microsatellite markers and phenotyping (growth and carbon uptake-related traits).

Clusters	Npure	Nww	Ncw	Ndom	Nno cluster	Total measured for phenotypic traits	Population name
Q1 (light green)	32	46	7	0	10	92	FR-W
Q2 (yellow)	0	52	5	0	4	57	FR-E
Q3 (pink)	28	1	0	0		28	FR-Lor
Q4 (blue)	61	21	1	0	6	85	DA
Q5 (purple)	23	4	3	0	4	34	IT
Q6 (dark green)	66	17	1	0	1	83	AUT
Q7 (red)	77	34	5	0	0	113	RO
Q8 (black; M. domestica)	40	0	46	21	8	73	DOM
Total	287	175	68	21	33	551 (584)
Total measured for height and number of leaves	282	167	63	21	32	533 (565)
Total measured for leaf chlorophyll and flavonol contents, NBI	129	82	22	6	18	239 (257)

Abbreviations: NBI, nitrogen balance index; N_pure, the number of seedlings assigned to a wild gene pool with a membership coefficient > 0.9; N_ww, the number of wild–wild hybrids (i.e. seedlings with a membership coefficient > 0.1 to a wild gene pool other than its own wild gene pool and a membership coefficient < 0.1 to the Malus domestica gene pool); N_cw, the number of crop–wild hybrids (i.e. seedlings assigned to the M. domestica gene pool with a membership coefficient > 0.1); N_{no cluster}, seedlings that could not be assigned to any defined gene pool. The total measured for phenotypic traits is the number of individuals measured for each phenotypic trait and included in the statistical analyses; the number in parentheses represents the initial sample size before data were filtered for statistical analyses. ‘Wild population name’ specifies populations defined with STRUCTURE for K = 8 excluding crop–wild hybrids and seedlings from misidentified mother trees (i.e. including only wild pure and wild–wild hybrids). Populations are as follows: Austria (AUT), Denmark (DA), Eastern France (FR-E), Lorraine in France (FR-Lor), Western France (FR-W), Italy (IT), Romania (RO) and M. domestica (DOM).

In the last week of the experiment, the superficial flavonol and chlorophyll contents were measured, and the nitrogen balance index (NBI) was calculated for three leaves per seedling. The superficial chlorophyll content is the chlorophyll concentration in the leaf epidermis (in micrograms per square centimetre), and the superficial flavonol content is an index of the flavonoid concentration (in micrograms per square centimetre) in this upper layer and is related to phenol accumulation and ultraviolet protection. The NBI assesses the nitrogen status of the leaf by calculating the ratio of chlorophyll to flavonols (related to nitrogen/carbon allocation) and is currently used as a proxy to estimate the foliar nitrogen content (Demotes-Mainard et al., 2008). Leaf chlorophyll, flavonol content and NBI are parameters correlated with plant carbon uptake via photosynthesis. Flavonol is a phenolic compound that contributes to plant vigour, acclimatization and adaptation to environmental constraints through various mechanisms, including its antioxidant activity. These traits were measured using a portable Dualex^® device (Force-A, Orsay, France), which uses a combination of fluorescence signals at various excitation bands to quantify pigments and chemical compounds and has already been calibrated on the apple tree (Hamann et al., 2018). Given that carbon uptake-related traits must be measured on the same day, a subsample of 257 of the 565 seedlings (Table 1; numbers in parentheses) was measured because of time limitations. Seedlings measured for carbon uptake-related traits were selected such that at least one seedling per mother tree and three seedlings per geographical site were sampled.

DNA extraction, microsatellite genotyping and genetic ancestry of the seedlings

At the end of the experiment, the leaves of each seedling were sampled for microsatellite genotyping. Genomic DNA was extracted using the NucleoSpin Plant DNA Extraction Kit II (Macherey & Nagel, Düren, Germany), according to the manufacturer’s instructions. Microsatellites were amplified by multiplex PCR with the Multiplex PCR Kit (QIAGEN, Inc.). We used 13 microsatellite markers (Ch01f02, Ch01f03, Ch01h01, Ch01h10, Ch02c06, Ch02c09, Ch02c11, Ch02d08, Ch03d07, Ch04c07, Ch05f06, GD12 and Hi02c07) in four multiplexes [MP01–MP04 (Cornille et al., 2012)].

Polymerase chain reaction was performed in a final reaction volume of 15 μL (7.5 μL of QIAGEN Multiplex Master Mix, 10–20 μm of each primer, with the forward primer labelled with a fluorescent dye, and 10 ng of template DNA). We used a touch-down PCR program (initial annealing temperature of 60 °C, decreasing by 1 °C per cycle to 55 °C). Genotyping was performed at the GENTYANE platform (INRAE Clermont-Ferrand) on an ABI PRISM X3730XL, with 2 mL of GS500LIZ size standard (Applied Biosystems). Alleles were scored with the GENEMAPPER v.4.0 software (Applied Biosystems). We retained only multilocus genotypes presenting <10 % of missing data.

Clones or closely related individuals can bias inferences about the population structure. We estimated the kinship coefficient between pairs of individuals (F_ij) with SPAGeDI 1.5d (Loiselle et al., 1995; Hardy and Vekemans, 2002) and removed individuals with high genetic relatedness with F_ij > 0.5.

The individual-based Bayesian clustering method implemented in STRUCTURE v.2.3.3 (Pritchard et al., 2000) was used to estimate the admixture between M. domestica and M. sylvestris, in addition to the population genetic structure of M. sylvestris. STRUCTURE uses Markov chain Monte Carlo simulations to infer the proportion of ancestry of genotypes from K distinct clusters. The underlying algorithm minimizes deviations from Hardy–Weinberg, and linkage disequilibria. The value of K ranged from one to ten. Ten independent runs were carried out for each value of K, and 500 000 Markov chain Monte Carlo (iterations after a burn-in of 50 000 steps were used. Clumpak (greedy algorithm) (Kopelman et al., 2015) was used to identify distinct modes in the ten replicated runs for each value of K. STRUCTURE analyses were run for the entire dataset (n = 584), plus 40 M. domestica genotypes from Western Europe (Supplementary data Table S2) included as a reference for the cultivated apple gene pool (Cornille et al., 2013b). The R package pophelper v.2.3.0 was used (Francis, 2016) to visualize bar plots. The amount of additional information explained by increasing K was determined using the ΔK statistic (Evanno et al., 2005), as implemented in STRUCTURE HARVESTER (Earl and vonHoldt, 2012). However, ΔK provides statistical support for the strongest but not the finest population structure (Puechmaille, 2016). Natural populations can display a hierarchical genetic structure, with a fine-scale population structure. Visual inspection of the bar plots was used to identify the value of K for which all clusters have well-assigned individuals and where additional clusters at higher values of K do not have well-assigned individuals. Therefore, the K value was the finest one, which can be higher than the K value of the strongest population structure identified by ΔK.

We used the membership coefficients at the best K value (inferred with STRUCTURE) for the crop and wild seedlings and defined the genetic status of each seedling. Initially, we separated cultivated apples (seedlings from mother trees that were misidentified in the field, i.e. with a membership proportion to the M. domestica gene pool, P_dom, > 0.9, referred to as ‘dom’ hereafter), crop–wild hybrids (seedlings with 0.1 > P_dom > 0.9, referred to as ‘cw’ hereafter) and pure M. sylvestris (seedlings with a membership coefficient to a given wild apple cluster > 0.9 and with P_dom < 0.1, referred to as ‘pure’ hereafter). We then separated ‘pure’ wild seedlings (i.e. seedlings with membership coefficients to a given wild gene pool > 0.9) from wild–wild admixed seedlings (i.e. seedlings with a membership coefficient to a given wild apple cluster < 0.9, referred to as ‘ww’ hereafter). Two effects were then tested using the statistical models described below: the genetic status effect (i.e. dom, cw, ww and pure) and the wild apple population effect (i.e. corresponding to the ‘pure’ populations detected with STRUCTURE). A seedling that could not be assigned to any cluster (with a membership coefficient to any cluster < 0.5) was removed from further analysis.

We computed descriptive population genetic estimates for each population (i.e. each cluster inferred with STRUCTURE, excluding admixed individuals with a membership coefficient < 0.9). Heterozygosity (expected and observed), Weir and Cockerham F-statistics and deviations from Hardy–Weinberg equilibrium were calculated with genepop v.4.2 (Raymond and Rousset, 1995; Rousset, 2008).

Phenotypic traits and correlations

Height, leaf growth rate and carbon-related traits were estimated for each seedling. The absolute height and leaf growth rates [AGR (Radford, 1967)], relative growth rates [RGR (Briggs et al., 1920)] and whole AGRs were estimated as follows:

$${\displaystyle \begin{array}{l}{\displaystyle \mathrm{A}\mathrm{G}\mathrm{R}\left(\mathrm{c}\mathrm{m}\mathrm{d}\mathrm{a}{\mathrm{y}}^{-1}\right)=\frac{(\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}{\mathrm{t}}_{t+1}-\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}{\mathrm{t}}_t)}{(\mathrm{d}\mathrm{a}\mathrm{t}{\mathrm{e}}_{t+1}-\mathrm{d}\mathrm{a}\mathrm{t}{\mathrm{e}}_t)}}\end{array}}$$

(1)

$${\displaystyle \begin{array}{l}{\displaystyle \mathrm{R}\mathrm{G}\mathrm{R}\left(\mathrm{c}\mathrm{m}\mathrm{d}\mathrm{a}{\mathrm{y}}^{-1}\mathrm{d}\mathrm{a}{\mathrm{y}}^{-1}\right)=\mathrm{A}\mathrm{G}\mathrm{R}/\mathrm{d}\mathrm{a}\mathrm{t}{\mathrm{e}}_t}\end{array}}$$

(2)

$${\displaystyle \begin{array}{l}{\displaystyle \mathrm{W}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{A}\mathrm{G}\mathrm{R}\left(\mathrm{c}\mathrm{m}\mathrm{d}\mathrm{a}{\mathrm{y}}^{-1}\right)=\frac{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}{\mathrm{t}}_{\mathrm{e}\mathrm{n}\mathrm{d}}-\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}{\mathrm{t}}_{\mathrm{b}\mathrm{e}\mathrm{g}\mathrm{i}\mathrm{n}\mathrm{n}\mathrm{i}\mathrm{n}\mathrm{g}}}{\mathrm{d}\mathrm{a}\mathrm{t}{\mathrm{e}}_{\mathrm{e}\mathrm{n}\mathrm{d}}-\mathrm{d}\mathrm{a}\mathrm{t}{\mathrm{e}}_{\mathrm{b}\mathrm{e}\mathrm{g}\mathrm{i}\mathrm{n}\mathrm{n}\mathrm{i}\mathrm{n}\mathrm{g}}}}\end{array}}$$

(3)

Note that for whole AGRs, the beginning of the experiment corresponded to days 7 and 11 for leaf and height measurements, respectively, and the last measurement was performed on day 60. The internode ratio, which represents the ratio between the number of leaves and the height of the seedling at day 60, was also considered as a fitness trait, because this value plays an essential role in apple tree architecture (Ripetti et al., 2008).

Seven phenotypic traits were therefore calculated for the entire dataset (565 seedlings; Table 1): height_AGR, height_RGR, whole_height_AGR, leaf_AGR, leaf_RGR, whole_leaf_AGR and internode. In addition, chlorophyll (Chl) content, flavonol (Flav) content and NBI were measured in a subsample of 257 seedlings (Table 1).

A preliminary exploration of the variation among phenotypic traits was performed. We used the cor R function to assess the correlation between phenotypic traits and seed sizes (length, width and weight) of the mother trees, and between phenotypic traits. The corrplot R package was used to visualize the correlations (Wei and Simko, 2017). Principal component analyses of the phenotypic traits of the seedlings were conducted with the FactoMineR R package (Lê et al., 2008).

Statistical analyses of phenotypic variation

A previous study demonstrated that crop-to-wild gene flow impacts early-stage growth rate (Feurtey et al., 2017). The effect of the genetic status of seedlings (dom, cw, ww or pure; Table 1) on fitness variation among seedlings was therefore tested. A linear mixed model was fitted to the data as follows:

$${\displaystyle \begin{array}{l}{\displaystyle \begin{array}{lll}& & {\displaystyle Y_{ijkl}\sim \mu +\mathrm{w}\mathrm{i}\mathrm{l}\mathrm{d}\mathrm{p}\mathrm{o}\mathrm{p}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{i}{\mathrm{n}}_i+\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{u}{\mathrm{s}}_j}\\ & & {\displaystyle +\mathrm{w}\mathrm{i}\mathrm{l}\mathrm{d}\mathrm{p}\mathrm{o}\mathrm{p}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{i}{\mathrm{n}}_i*\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{u}{\mathrm{s}}_j}\\ & & {\displaystyle +\mathrm{m}\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}{\mathrm{r}}_k+{\mathrm{e}}_{ijkl}}\end{array}}\end{array}}$$

(4)

where Y_ijk is the phenotypic trait (height or leaf growth rate or a carbon uptake-related trait) of a seedling from mother k from population i with the genetic status j, wild population of origin is the fixed effect of the population of origin i of the seedling inferred with STRUCTURE, genetic status is the fixed effect of the genetic status of the seedling (i.e. pure, dom, ww or cw), the interaction between the two fixed effects, mother tree is a normally distributed random effect with its own mean and variance parameters and e_ijkl is the residual. The mother tree of each seedling was used as a random factor to avoid pseudo-replication owing to the presence of multiple half-siblings (i.e. from the same mother tree). We ran the model (4), but replaced the genetic status effect with the P_dom fixed effect. We gradually removed interactions and effects depending on their significance. In addition, we evaluated the differences in the effect on trait variation using a contrast analysis. We fitted the data to the model using the lme4 R package (Bates et al., 2015). The statistical significance of an effect was assessed with a type II Wald χ² test.

For phenotypic traits estimated from the number of leaves (leaf_AGR, leaf_RGR and whole_leaf_AGR), a log-link function was used, and the residual distribution was fitted to a negative binomial distribution (function glm.nb in R package lme4). For phenotypic traits defined from the height of the seedling (height_AGR, height_RGR and whole_height_AGR) and for chlorophyll content, flavonol content and NBI, a similar linear mixed model was run, but with a residual term that was assumed to be normally distributed.

Heritability estimates

The level of heritability estimate provides an indication of the potential extent of the response to the selection of a trait. Therefore, computing heritability can allow us to determine whether the traits measured in this study could be good candidates to respond to selection. Heritability estimates were calculated using only pure and wild–wild M. sylvestris seedlings (Table 1). We fitted each fitness proxy with a linear mixed model as follows:

$${\displaystyle \begin{array}{l}{\displaystyle Y_{ijk}=\mu +F_i+C_j+{\mathrm{e}}_{ijk}}\end{array}}$$

(5)

where Y_ijk is the fitness proxy (growth rate or carbon uptake-related trait) of the kth seedling from mother tree i, member of the jth genetic cluster, μ the overall fixed mean of the population, F_i the random effect of the ith mother tree, C_j is the fixed effect of the jth genetic cluster and e_ijk the random error term. The model was fitted using REML (restricted maximum likelihood). Calculations were performed by the lme function of the R library nlme (Pinheiro et al., 2018). The output of lme provides estimates for the variance components, the corresponding s.d., and the best unbiased linear predictors for random effects. Genetic parameters were then calculated as follows:

1. the additive genetic variance: $\mathrm{V}\mathrm{A}=4{\sigma }_{\mathrm{F}}^2$⁠, with ${\sigma }_{\mathrm{F}}^2$ representing the between-mother tree variance;

2. the corresponding coefficient of variation: $\mathrm{C}\mathrm{V}\mathrm{A}=\frac{\sqrt{\mathrm{V}\mathrm{A}}}{\mu }$⁠;

3. the phenotypic variation: $\mathrm{V}\mathrm{P}={\sigma }_{\mathrm{F}}^2+{\sigma }_{\mathrm{E}}^2$⁠, with ${\sigma }_{\mathrm{E}}^2$ representing the residual variance;

4. the corresponding coefficient of variation: $\mathrm{C}\mathrm{V}\mathrm{P}=\frac{\sqrt{\mathrm{V}\mathrm{P}}}{\mu }$⁠;

5. Narrow-sense heritability: $h^2=\mathrm{V}\mathrm{A}/\mathrm{V}\mathrm{P}$⁠; and

6. Dickerson’s approximation of s.d: $\mathrm{s}.\mathrm{d}.\left(h^2\right)\approx \frac{4\mathrm{s}.\mathrm{d}.\left({\sigma }_{\mathrm{F}}^2\right)}{\mathrm{V}\mathrm{P}}$⁠.

Test for isolation by ecology

Only pure and wild–wild hybrid M. sylvestris seedlings were sampled for the IBE analysis (n = 449, 21 sites; Table 1). The IBE pattern, i.e. the contribution of climatic distance (isolation by climate) and phenotypic distance [isolation by adaptation (IBA)] to the genetic structure, taking into account geographical distance, was evaluated using distance-based redundancy analysis (db-RDA), which can be used when the response variable is a distance matrix [here, it is a genetic distance matrix (F_ST) across 21 sampled sites] and the explanatory variables are in vector form as follows. The geographical distance between sampled sites, estimated with SPAGeDI 1.5d (Hardy and Vekemans, 2002), which underlies an IBD process, was represented by vectors with positive eigenvalues of a principal coordinate of a neighbour matrix (PCNM) (Borcard and Legendre, 2002). A total of 19 bioclimatic variables, downloaded from the Worldclim2 database (30″ resolution; https://www.worldclim.org/data/worldclim21.html) and representing annual and seasonal trends and extremes averaged over the years 1970–2000 and averaged for the Pleistocene period (20 000 years ago) (Gamisch, 2019), were used to test for an isolation-by-climate pattern. Growth rate (n = 551; Table 1) and carbon uptake-related traits (n = 239; Table 1) averaged per site were used to test for an IBA pattern.

To identify the variables that explained the genetic structure of M. sylvestris, a db-RDA using the ‘capscale’ function (Oksanen et al., 2014) was run on a complete model that included all investigated variables (i.e. PCNM components, growth rates, carbon uptake-related traits and 19 bioclimatic variables). The best variables were selected for an optimum model with the function ‘step’ based on the Akaike information criterion (AIC). Given that db-RDA does not provide information on the relative contribution of each variable of the model, a variance partitioning analysis was run using the ‘varpart’ function from the R package ‘vegan’ (Peres-Neto et al., 2006).

RESULTS

Genetic status of seedlings

No clones or closely related individuals were detected (Supplementary data Fig. S1); therefore, 584 seedlings were included in the STRUCTURE analyses.

STRUCTURE revealed a clear spatial population genetic structure of M. sylvestris in Europe and crop–wild admixture (Fig. 1; Supplementary data Fig. S2). From K = 2 to K = 8, several clusters appeared. When K > 8, STRUCTURE did not reveal any further substructures but only additional clusters with highly admixed individuals (Supplementary data Fig. S2). Therefore, although the ΔK indicated that the most likely value of K was five (Supplementary data Fig. S3), K = 8 was the finest population structure and was retained in subsequent analyses. For K = 8, we found that the M. domestica reference varieties were admixed with Italian and Western French M. sylvestris. Note that cultivars used as references originated from Western Europe (Supplementary data Table S2). Conversely, we detected 68 M. sylvestris seedlings with P_dom > 0.1 (considered to be cw hybrids), corresponding to 11.6 % of the seedlings (n = 584; Fig. 1; Table 1; Supplementary data Figs S4–S6). We also found 21 seedlings with a membership coefficient to the M. domestica gene pool > 0.9, corresponding to 4 % of the seedlings. Nearly all Spanish seedlings were assigned to the M. domestica gene pool with membership coefficients > 0.1 (i.e. 26 cw hybrids and 13 individuals assigned to the M. domestica gene pool) and showed admixture only with the wild Italian purple gene pool (Supplementary data Figs S4 and S6). A total of 33 individuals could not be assigned to any cluster (i.e. individuals with a membership coefficient to any cluster < 0.5).

Fig. 1.

Bayesian clustering of Malus sylvestris seedlings sampled in this study (n = 584) and reference samples of Malus domestica (n = 40) inferred with STRUCTURE at K = 8 and its associated map of mean membership per sampled site. Each individual is represented by a vertical bar partitioned into clusters. Visualization was improved by sorting genotypes by country. Countries are separated by a white line. The M. domestica reference samples are shown on the far left of the map in the Atlantic. Circle size is proportional to the number of individuals within the cluster (scale shown in the top right-hand corner).

We, therefore, identified 68 cw, 21 dom, 167 ww and 282 pure seedlings (Table 1; n = 551). After removal of cw hybrids (n = 68), seedlings sampled from misidentified mother trees (n = 21), the M. domestica reference samples (n = 40) and individuals with a membership coefficient to any cluster < 0.5, the following seven wild apple populations (i.e. groups of seedlings with a membership coefficient > 0.5 to a wild apple cluster) were defined: French Western (FR-W; n = 77), French Eastern (FR-E; n = 50), French Lorraine (FR-Lor; n = 28), Danish (DA; n = 78), Italian (IT; n = 27), Austrian (AUT; n = 81) and Romanian (RO; n = 108) (Supplementary data Fig. S6). Each M. sylvestris population exhibited a high level of genetic variation (Supplementary data Table S3). The Romanian population was the most genetically differentiated and was close to the Austrian population; the Danish and French Western populations were the closest genetically (Supplementary data Fig. S7), which is congruent with previous results from 26 microsatellite markers (Cornille et al., 2015).

No effect of seedling genetic status on phenotypic variation

Variations and correlations among phenotypic traits are presented in the Supplementary data (Figs S8–S12). Significant correlations between seed length and weight of the mother tree and height of the seedling were observed (P < 0.01). Leaf growth rates (leaf_RGR, leaf_AGR and leaf_whole AGR) were significantly correlated, as were height growth rates (height_RGR, height_AGR and height_whole_AGR). Chlorophyll content was positively correlated with NBI, whereas flavonol content was negatively correlated with NBI. Heritability estimates were moderate to high for all traits except growth rates based on leaf number (AGR_leaf, RGR_leaf and whole_AGR_leaf; Supplementary data Table S4). Given the limited sample size, these estimates must be taken cautiously, as reflected by the large values of s.d.

We did not find any significant effect of the genetic status of the seedlings (i.e. pure, ww, cw or dom; Supplementary data Table S5; Fig. S13) or P_dom (Supplementary data Table S6; Fig. S14) on phenotypic traits, except for a marginally significant effect (P = 0.04; Supplementary data Table S6) of P_dom on the whole height AGR: seedlings with higher levels of introgression by M. domestica had higher whole height AGRs. We, therefore, removed the seedling genetic status and P_dom effects from the model (4), in addition to dom and cw individuals. Thus, the final model included only wild apple seedlings (i.e. pure and ww, n = 449) and focused on the wild population of origin effect (Table 2).

Table 2.

Final model depicting the effects of the Malus sylvestris population to which each seedling belonged (i.e. cluster inferred with STRUCTURE for K = 8) on phenotypic traits (i.e. height, number of leaves, internode, chlorophyll and flavonol contents and nitrogen balance index) measured on 533 individuals.

Explanatory variable	Population			Mother tree		Model			Mother tree
Fitness	χ2	P-value	d.f.	REML	s.d.	AIC	R2	Corrected R2	R2
Height_AGR	17.863	0.007***	6	1229	0.204	1248	0.047	0.091	0.044
Height_RGR	12.846	0.045*	6	−2264	0.006	−2245	0.041	0.147	0.106
Leaf_AGR	4.302	0.6359	6	134	0.036	152	0.01	0.028	0.018
Leaf_RGR	4.302	0.6359	6	134	0.036	152	0.01	0.028	0.018
Whole height AGR	22.243	1.00e-03***	6	630	0.175	650	0.074	0.192	0.118
Whole leaf AGR	36.326	2.38e-06***	6	−1277	0.009	−1258	0.09	0.118	0.028
Height	31.623	1.93e-05***	6	4113	11.69	4,131	0.119	0.301	0.182
Number of leaves	22.285	0.001***	6	–	0.084	2,659	0.052	0.064	0.012
Chlorophyll	14.418	0.025*	6	1181	1.352	1199	0.074	0.171	0.097
Flavonol	6.7752	0.342	6	−124	0.091	−105	0.043	0.299	0.256
NBI	1.838	0.934	6	1735	7.6617	1754	0.011	0.224	0.213
Internode (number of leaves/height)	17.768	0.007***	6	−1328	0.009	−1309	0.044	0.073	0.029

Abbreviations: AIC, Akaike information criterion; NBI, nitrogen balance index; –, model without any significant effect. ***P < 0.001; **0.01 < P < 0.001; *0.05 < P < 0.01.

Significant variation in growth rates and chlorophyll content among populations

The mean height variation during the experiment among seedlings from different populations is shown in Figure 2. There was significant variation among seedlings from other populations in certain growth-related traits (Table 2). On average, seedlings belonging to the Austrian population were taller (+11 cm, P = 0.047), whereas Romanian (−14.9 cm, P = 0.008) and Italian (−18.7 cm, P = 0.044) seedlings were shorter (Fig. 2; Supplementary data Fig. S15) than seedlings from other populations. Seedlings from other populations did not show a significant difference in height. In addition, the number of leaves and height traits were negatively correlated (r = −0.3, P < 0.001). The Austrian population presented the lowest number of leaves (average = 5, s.d. = 4), whereas seedlings belonging to the Romanian population had the highest number of leaves (average = 8, s.d. = 7; Supplementary data Figure S16). The Romanian population also had the largest internode (+0.02 leaves cm⁻¹, P = 0.024).

Fig. 2.

Cumulative height growth over time in apple seedlings in controlled conditions. Seedlings measured were pure and wild–wild hybrid Malus sylvestris seedlings (n = 449) and seedlings assigned to the Malus domestica gene pool (n = 21), as detected with STRUCTURE for K = 8. The 40 reference M. domestica individuals were not measured in controlled conditions and are therefore not shown here. Vertical lines represent the s.d. Populations: Austria (AUT, n = 81), Denmark (DA, n = 78), M. domestica (DOM, n = 21, which includes 13 Spanish genotypes and seedlings from other countries), Eastern France (FR-E, n = 50), Lorraine in France (FR-Lor, n = 28), Western France (FR-W, n = 77), Italy (IT, n = 27) and Romania (RO, n = 108).

Chlorophyll content differed among populations, with seedlings from the Italian population producing, on average, more chlorophyll (+4.14 µg cm⁻², P = 0.039; Supplementary data Fig. S17) than seedlings from other populations. Flavonol content and NBI did not differ significantly among populations.

Significant IBD and IBC

Correlation plots between bioclimatic variables are provided in Supplementary data Figs S18 and S19; however, all variables were included in the analysis, because db-RDA can cope with correlated variables. The optimal model was chosen according to its best AIC value. The optimal model explained ≤25.9 % of the genetic structure (adjusted R² = 69.9 %, P = 0.001) and contained seven variables (four geographical and three bioclimatic variables) (Table 3): the geographical distance is represented by the first, second, third and sixth axes of the PCNM analysis and three past climatic variables (Bio3, isothermality; Bio6, minimum temperature of the coldest month; and Bio 9, mean temperature of the driest quarter). In total, IBD explained 47 % of the variance of the wild apple tree population genetic structure, whereas IBC explained 22 % (Supplementary data Fig. S20). Taking geographical distance into account, we did not find a pattern of IBA, i.e. covariation between phenotype and genetic divergences.

Table 3.

Contribution of geography and climate to the genetic variation observed among Malus sylvestris seedlings. Distance-based redundancy analyses tested the effects of geography, climate and phenotype on the genetic differentiation (from 13 microsatellites) among 21 sites in the European crabapple. Only the significant variables are listed.

	db-RDA
	Percentage of variance explained	d.f.	P-value	Adjusted R2
Global analysis	25.9	7	0.001	69.9
Residuals	11.3	13	–
Geography (IBD, PCNM 1-2-3-6)	14.9	4	<0.015
Environment (IBC_LGM, BIO3_LGM, BIO6_LGM, BIO9_LGM)	11.04	3	<0.015
Residuals	11.3 %	–	–

Abbreviations: BIO3_LGM, isothermality (BIO2/BIO7) (×100); BIO6_LGM, minimum temperature of the coldest month; BIO9_LGM, mean temperature of the driest quarter; IBD, isolation by distance; IBC_LGM, isolation by climate during the Last Glacial Maximum.

DISCUSSION

This study assessed the relationship between phenotypic and genetic variation among populations of a wild contributor to the cultivated apple genome (Cornille et al., 2012), M. sylvestris, a CWR species that occurs naturally along a climatic gradient in Europe. Bayesian clustering revealed a substantial number of seedlings introgressed by M. domestica. With the hybrids removed, seven populations of M. sylvestris distributed across Europe were detected and showed phenotypic variation in growth and chlorophyll content. Based on the IBA pattern estimated from the phenotypic traits measured in this study, this phenotypic variation was not adaptive. However, the significant association between population genetic variation and the LGM climate suggests that the European crabapple might be locally adapted to the past climate conditions of the LGM. The results of this study, therefore, indicate the occurrence of adaptive divergence related to climate in a wild contributor to the cultivated apple genome. This might help us to use its diversity better, providing options for mitigating the impact of climate change on the cultivated apple through breeding (Warschefsky et al., 2014; Prohens et al., 2017; Satori et al., 2022).

Ongoing crop-to-wild gene flow in the European crabapple

We showed substantial gene flow from M. domestica to the European crabapple gene pool, with 11.6 % of seedlings, mainly from Western Europe, introgressed by M. domestica. Introgression rates were lower compared with previous studies [i.e. 37 % in the study by Cornille et al. (2013b) and 23.1 % in the study by Cornille et al. (2015)]. However, these studies genotyped more mother trees (i.e. n = 756 and n = 1889, respectively), which could explain the difference in estimates of crop-to-wild gene flow. Crop-to-wild gene flow is therefore still ongoing in the European crabapple. The higher number of seedlings from Western Europe that are introgressed by M. domestica can be explained by the use of reference Western European cultivated apple varieties. Alleles specific to Eastern and Northern cultivated apple varieties could be missed and might decrease the probability of detecting crop-to-wild introgression events in non-Western wild populations (from Eastern and Northern Europe). Note that the Spanish seedlings sampled in this study were the progeny of trees growing in a location known to have high levels of introgression by M. domestica (G. Alins, IRTA, Lleda, Spain, pers. comm.). It is even possible that the mother trees of these seedlings were M. domestica and not M. sylvestris.

The consequences of crop-to-wild introgression on phenotypic variation between crop and wild individuals are poorly understood in perennial fruit trees. One study has shown that crop–wild hybrid apple seedlings have higher growth rates and germinate earlier than wild apple seedlings (Feurtey et al., 2017). We did not detect any effect of the genetic status of a seedling (pure, ww, cw or dom) or the level of introgression (P_dom) on growth and carbon uptake-related fitness proxies. This could be attributable to the low number of samples from the cw and dom categories.

Signs of local adaptation to past climate in the European crabapple

In controlled conditions, seedlings from different populations exhibited significantly different growth and physiology. Seedlings belonging to the Austrian population were the tallest, with the highest absolute growth rate and the lowest number of leaves; in contrast, Romanian seedlings were the shortest and had the lowest absolute growth rate and the highest number of leaves. Italian seedlings had the highest chlorophyll content. We therefore tested whether this phenotypic variation was adaptive. For instance, the seedlings from the Austrian population might be fitter in the climatic conditions simulated in this controlled condition experiment. However, considering the geographical distance, we found no significant covariation between genetic and phenotypic variation. This suggests a lack of divergent selection on carbon uptake or growth traits, which are often associated with plant responses to climate (Bussotti et al., 2015). Therefore, the phenotypic variation we observed among populations in controlled conditions might result from genetic drift alone and not from divergent selection. Alternatively, although the specific traits we selected are among the traits that are generally considered to be related to responses to climate (Kühn et al., 2021), they might not be perfect candidates for investigating divergent selection by climate. Another explanation for the lack of observation of divergent selection on carbon uptake or growth traits could be that we did not phenotype enough seedlings from each genetic group. We observed a high variation in each phenotypic trait and their heritability estimates, suggesting that the traits we studied might be relevant but that a larger number of seedlings could be phenotyped and analysed. However, some studies have found that even with large sample sizes, the s.e. of heritability estimates can still be large and vary significantly between experimental designs (Visscher and Goddard, 2015). The reasonably high heritability estimates for most traits considered here could be consistent with relatively weak within-population selection, enabling the maintenance of ample additive genetic variation (Wheelwright et al., 2014). Furthermore, high variation in seedling traits combined with high heritability estimates could suggest a substantial amount of genetic variation for adaptation to work on. Altogether, the effect of the population of origin on phenotypic traits suggests a genetic basis for this variation. However, we need to increase the sample size and measure new traits to draw more precise conclusions on the occurrence of adaptive phenotypic variation in the European crabapple.

Although the association between phenotypic and genetic variation was not significant we observed a significant association between LGM climate and genetic variation in the European crabapple, suggesting local adaptation. Given that M. sylvestris is distributed across various climatic conditions, we indeed investigated the role of climate in shaping the genetic variation among populations of the European crabapple without considering phenotypic trait variation. We tested for an IBE pattern, where the pattern of neutral genetic variation covaries with ecological variables (here, climate). There was no combined effect of geographical and climatic distance (IBD ∩ IBC), which allowed us to assess the contribution of these processes separately (Wang and Bradburd, 2014). We showed that IBD and IBC played a significant role (adjusted R² = 47 % and adjusted R² = 22 %, respectively) in the genetic differentiation of European crabapple populations. Weak but significant IBD has been identified previously in wild apple relatives of the cultivated apple (i.e. M. sylvestris, M. orientalis and M. sieversii; Cornille et al., 2013a, b; Cornille et al., 2015). Here, we used 13 of the 26 microsatellite markers used in previous studies; the lack of resolution of the 13 simple sequence repeats might explain the lack of IBD. However, this is unlikely, because the Bayesian inferences method was previously able to detect genetic clusters in wild apples (Cornille et al., 2013a, b; Cornille et al., 2015). A weak IBD pattern suggests that M. sylvestris has high dispersal capacities (Coart et al., 2006; Larsen et al., 2006; Cornille et al., 2013a; Cornille et al., 2015; Feurtey et al., 2017; Reim et al., 2017). A weak IBD is explained by a self-incompatibility system that prevents self-fertilization (Brown, 1992), pollen dispersal by bees, beetles and flies, and endozoochorous seed dispersal by large mammals, such as ungulates, wild boars, brown bears or humans (Larsen et al., 2006).

We showed that in addition to IBD, IBC persisted after considering the geographical distance. Climate can impose divergent selection pressures on different locations and thus reduce gene flow between populations. For instance, divergent selection imposed by climate can limit the reproductive success of individuals moving between different climates for which they are adapted, such that IBC contributes to genetic differentiation among populations (Wang and Bradburd, 2014). The main variables explaining genetic differentiation in the European crabapple were related to temperature during the LGM. This suggests that the European crabapple might be adapted to its past climate but not to its current one. In wind-dispersed trees, local adaptation to the current climate has been demonstrated (Savolainen et al., 2013; Kremer and Hipp, 2019; Pyhäjärvi et al., 2020), but to our knowledge, no study has shown local adaptation to past climate conditions in an insect-pollinated tree species.

Factors other than climate can also shape the adaptive divergence between populations. Local adaptation to biotic factors, such as the presence of other species, is possible. Malus sylvestris is a species that needs high levels of light and is not very competitive. Some of the observed variations of the seedlings in growth rate, number of leaves and chlorophyll content might derive from adaptations that would be advantageous within the local niche. In Romania, the populations sampled were in forest edges, middle succession woodlands (with hawthorn and wild pear) or grasslands; no wild apples were found in mature forests. There, M. sylvestris does not compete with other woody species but with grasses and shrubs in the seedling phase; this might explain the slower growth in height and the larger number of leaves (leading to a larger leaf area and more shading of the competing grasses). In France, the trees are present at the edge of mature forests, and in Austria, in mature forests. It is hard to draw clear conclusions; further investigations of the ecology of the wild populations in situ are needed. The rhizosphere composition among seedlings from different populations can impact the phenotypic variation and, potentially, plant fitness in response to climate (Trivedi et al., 2022). Here, seeds were cleaned with chlorine to avoid the effect of local micro-organism community and seed growth. Still, investigations of phenotypic variation among cleaned and uncleaned seeds can also help to assess the role of the rhizosphere in the divergence of the wild apple population. Local adaptation of fruit trees to biotic factors, including parasites (Olvera-Vazquez et al., 2021), also deserves further investigation. Besides selection, the role of phenotypic plasticity in enabling growth and optimal fitness in changing environments also needs to be evaluated carefully (Benito Garzón et al., 2011).

Further investigations are needed on local adaptation and phenotypic plasticity in response to climate in the European crabapple

Our study raises questions regarding the response of wild apple populations to climate change. More data are needed to draw clear conclusions. The adaptation of tree species to climate is complex (Bussotti et al., 2015). For instance, in one study, Eucalyptus camaldulensis, variation in leaf traits and performance was unrelated to the climate of genotype provenance (Asao et al., 2020), whereas in another part of Australia, the same species displays variation in several photosynthetic traits that are related to the climate of genotype provenance (Dillon et al., 2018). In contrast, collective differences in leaf morphology and photosynthetic physiology associated with the length of the growing season, temperature and the level of insolation in several Populus species might be adaptive (Keller et al., 2011; Kaluthota et al., 2015). Further investigations on local adaptation and phenotypic plasticity in response to climate in the European crabapple are needed. Genomic data will help in determining the relative influence of adaptive and neutral processes on climate-driven divergence by scanning the genomes of trees from different populations in Europe. Comparing the fitness of seedlings from different populations in reciprocal transplants (in controlled or natural conditions) will also allow us to investigate local adaptation and phenotypic plasticity. Further studies using additional genetic markers (single nucleotide polymorphism) and measuring phenotypic traits in different climatic conditions are therefore needed. This study is, nevertheless, a starting point for future breeding and conservation programmes of a CWR of an emblematic temperate fruit tree, because it characterized the phenotypic and genetic variation of seedlings in the wild that can be used as ex situ sources to enrich the crop gene pool. Some of the seedlings included have sequenced genomes and have been planted in several orchards in France (https://www.ideev.universite-paris-saclay.fr/en/the-orchard) and are measured each year for several phenotypic traits. Information and samples from these orchards can be requested from the corresponding author. Genetic variation of trees in these ex situ orchards can help to enrich the cultivated apple gene pool. Certain traits related to adaptation to climate change can be introduced into future apple varieties (Warschefsky et al., 2014; Prohens et al., 2017; Satori et al., 2022).

SUPPLEMENTARY DATA

Supplementary data are available at Annals of Botany online and consist of the following.

Table S1: details of Malus sylvestris material used in this study. Table S2: Malus domestica reference samples used in this study. Table S3: genetic diversity estimates for the seven populations of Malus sylvestris inferred at K = 8 with STRUCTURE. Table S4: variance and heritability of the traits measured in 533 individuals. Table S5: effects of status of the apple seedling and of the population to which they belong on phenotypic traits measured on 533 individuals. Table S6: effects of admixture proportion of a wild genotype assigned to the Malus domestica gene pool and of the population to which they belong on phenotypic traits measured on 533 individuals. Figure S1: distribution of kinship between pairs of Malus sylvestris and Malus domestica individuals. Figure S2: Bayesian clustering for Malus sylvestris and Malus domestica, inferred from K = 2 to K = 9 with STRUCTURE v.2.3.3. Figure S3: ΔK estimated for each K cluster inferred with STRUCTURE v.2.3.3 for Malus sylvestris and Malus domestica. Figure S4: distribution of crop–wild hybrids for the apple tree in Europe inferred with STRUCTURE v.2.3.3 at K = 8. Figure S5: distribution of the proportion of admixture of seedlings used in this experiment with the Malus domestica gene pool, P_dom, inferred with STRUCTURE at K = 8. Figure S6: distribution of the genetic status of the seedlings inferred at K = 8 with STRUCTURE, summed up per geographical site. Figure S7: pairwise genetic differentiation between populations of pure wild Malus sylvestris seedlings inferred at K = 8 with STRUCTURE. Figure S8: Pearson correlations between seed size of the mother trees and plant height and the number of leaves of the seedlings in Malus sylvestris. Figure S9: principal component analysis representing variation among phenotypic traits. Figure S10: correlation plot among phenotypic traits measured on the apple seedlings in controlled conditions. Figure S11: principal component analysis among calculated growth rate variables. Figure S12: correlation plot among growth rates obtained using cor() R-function for seedlings used in this study. Figure S13: evolution of mean height of seedling over time in growth chamber depending on the status defined with STRUCTURE for K = 8. Figure S14: the mean absolute growth rate of seedlings plotted against their level of Malus domestica ancestry. Figure S15: variation in height among the eight populations detected with STRUCTURE at K = 8 for seedlings of Malus sylvestris grown in controlled conditions. Figure S16: variation in the number of leaves among the eight populations detected with STRUCTURE at K = 8 for seedlings of Malus sylvestris grown in controlled conditions. Figure S17: variation in chlorophyll content among the eight populations detected with STRUCTURE at K = 8 for the seedling of Malus sylvestris grown in controlled conditions. Figure S18: correlation plot among the 19 bioclimatic variables averaged over 1970–2000 used to test for isolation by climate in Malus sylvestris. Figure S19: correlation plot among the 19 bioclimatic variables averaged for the Pleistocene period used to test for isolation by climate in Malus sylvestris. Figure S20: variance partitioning analysis of db-RDA results obtained for Malus sylvestris.

ACKNOWLEDGEMENTS

We are grateful to the INRAE MIGALE bioinformatics platform (http://migale.jouy.inra.fr) for providing help and support, particularly Véronique Martin, Eric Montaubon and Valentin Loux. We also thank Adrien Falce, Olivier Langella and Benoit Johannet for help and support on the INRAE Génétique Quantitative et Evolution – Le Moulon lab cluster and the Plateforme de Génotypage GENTYANE INRAE UMR 1095. We thank Einar Flensted-Jensen, Michel Bonnessee, Alberto Dominici, Jean-Luc Ricaud, Eric Wolff, Alain Guinet and Line Caminade, Richard Le Menn and Bernard Gambier for providing seeds. USDA is an equal opportunity provider, employer and lender. Mention of trade names or commercial products is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the US Department of Agriculture. A.C., G.M.V., A.R., S.B. and T.M.U. conceived and designed the experiments; A.C., G.M.V., A.F., T.M.U. and T.K. obtained the funding; A.C., G.M.V., A.F., T.M.U., K.A., S.O.-V., R.R., X.C., T.K. and C.R. sampled the material; X.C., C.R., A.V., A.R., G.L., K.A.O., R.R., M.L.G., H.B., V.C., H.C., S.O.-V. and M.F. performed the molecular biology analyses; A.C., A.F., K.A., and X.C. analysed the data. The manuscript was written by A.C., K.A., and A.F., with essential input from other co-authors.

FUNDING

A.R. and T.M.U. were funded, in part. through project PN2019-2022/19270201—BIODIVERS 3. We thank the Laboratoire d’Excellence BASC (Biodiversité, Agroécosystèmes, Société, Climat), the Institut Diversité Écologie et Évolution du Vivant (IDEEV), the European LEADER program and ATIP-Avenir CNRS Inserm for funding.