# Investigating the genetic architecture of biotic stress response in stone fruit tree orchards under natural infections with a multi-environment GWAS approach

Marie Serrie, Vincent Segura, Alain Blanc, Laurent Brun, Naïma Dlalah, Frédéric Gilles, Laure Heurtevin, Mathilde Le Pans, Véronique Signoret, Sabrina Viret, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marie Serrie, Vincent Segura, Alain Blanc, Laurent Brun, Naïma Dlalah, et al. Investigating the genetic architecture of biotic stress response in stone fruit tree orchards under natural infections with a multi-environment GWAS approach. 2024. hal-04759062

HAL Id: hal-04759062

https://hal.inrae.fr/hal-04759062v1

Preprint submitted on 29 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Investigating the genetic architecture of biotic stress response in stone fruit tree orchards under natural infections with a multi-environment GWAS approach

Marie Serrie¹, Vincent Segura²,³, Alain Blanc⁴, Laurent Brun⁵, Naïma Dlalah¹, Frédéric Gilles¹, Laure Heurtevin¹, Mathilde Le-Pans⁵, Véronique Signoret¹, Sabrina Viret¹, Jean-Marc Audergon¹, Bénédicte Quilot¹, Morgane Roth¹*

¹INRAE, UR GAFL, Avignon, France

²AGAP Institut, Université Montpellier, CIRAD, INRAE, Institut Agro, Montpellier, France

³Geno-Vigne®, IFV, INRAE, Institut Agro, Montpellier, France

⁴INRAE, UE A2M, Avignon, France

⁵INRAE, UERI Gotheron, Saint-Marcel-Lès-Valence, France

Morgane Roth

INRAE GAFL

67 Allée des Chênes, 84140 Avignon, France

email: morgane.roth@inrae.fr

**Running title:** Peach and apricot and biotic stress low susceptibility

<footnote>*Corresponding author:</footnote>
**Abstract**

The mapping and introduction of sustainable plant immunity to pests and diseases in fruit tree is still a major challenge in modern breeding. This study aims at deciphering the genetic architecture underlying resistance or tolerance across environments for major pests and diseases in peach (P. persica) and apricot (P. armeniaca). We set up a multi-environment trial (MET) approach by studying two core collections of 206 peach and 150 apricot accessions deployed under low phytosanitary conditions in respectively three and two environmentally contrasted locations in South-East of France. To capture the complex dynamics of pest and disease spread in naturally infected orchards, visual scoring of symptoms was repeated within and between 3 years, for five and two pests and diseases respectively for peach and apricot, resulting in the maximum of damage score and the AUDPC. These traits were used as phenotypic inputs in our genome-wide association studies (GWAS) strategy, and leading to the identification of: i) non-additive genotype-phenotype associations, ii) environment-shared QTLs iii) environment-specific QTLs, and iv) interactive QTLs which changes in direction ('antagonist') or intensity ('differential') according to the environment. By conducting GWAS with multiple methods, we successfully identified a total of 60 high confidence QTLs, leading to the identification of 87 candidate genes, the majority belonging to the Leucine-rich repeat containing receptors (LRR-CRs) family gene. Finally, we provided a comparative analysis of our results on peach and apricot, two closely related species. The present results contribute to the development of genomics-assisted breeding to improve biotic resilience in *Prunus* varieties.

**Key words**

Biotic stress; genome-wide association studies (GWAS); multi-environment trials (METs); Genotype × Environment interactions; resistance; tolerance; stone fruit tree; low phytosanitary inputs; core collections; *Prunus*
Introduction

Biotic stresses represent a major threat to agricultural production with a massive diversity of organisms being able to affect plants, including a large range of viruses, bacteria, fungi, phytoplasma, nematodes and pests [1,2]. Reducing the impact of these biotic stresses in a sustainable and environment-friendly way is urgent and, as such, reducing plant susceptibility through genetic progress appears as an effective and promising strategy [3]. It has been long recognized that host defence against parasites and pathogens can be divided into two conceptually different mechanisms: the ability to reduce and control pathogen (resistance) and the ability to limit the damage caused by a given within-host pathogen load without necessarily reducing this pathogen load (tolerance) [4,5]. Developing cultivars that present durable low susceptibility to biotic stresses is therefore an important breeding goal, in particular for perennial woody crops, given their long generation times and lifespans in orchards [6]. After multiple episodes of resistance breakdown due to the introgression into crop cultivars of monogenic resistances [7,8], research is now focusing on polygenic defence traits which are likely more durable but also more complex to understand and to introduce in elite lines [9]. The combination of major resistance genes with quantitative resistance/tolerance factors have already shown promising results as an alternative to ensure greater durability of crop protection [10,11]. In addition to specific host resistances genes, broad-spectrum defence factors against multiple diseases represent another promising solution for sustainable plant defence [12]. Some degree of tolerance has already been observed in different landraces cultivar and wild relatives, as *P. davidiana*, a relative of peach (*P. persica*) carrying resistance genes against green peach aphids, powdery mildew and *Plum pox virus* (PPV) [13–15]. In this context, the natural genetic diversity of crops is a valuable source of untapped favourable allelic combinations which deserves to be explored [16].

Genome Wide Association Studies (GWAS) have been extensively used to understand the genetic basis of complex traits such as quantitative disease resistance or tolerance [17]. This whole-genome approach exploits ancestral linkage disequilibrium among diverse germplasm to detect associations between phenotypic traits and polymorphisms, also called quantitative trait loci (QTL) [18,19]. QTLs identified by GWAS can be used in breeding programs through the development of marker-assisted selection (MAS) in relevant germplasm [20,21]. Yet, detecting QTLs related to disease resistance or tolerance can be challenging since phenotypic variation is often modulated by environmental effects and interactions between genotype and environment (G×E) [22–24]. The consequences of G×E can be illustrated by a change in the phenotypic ranking of individuals across environments which could be explained by a modulation of the size and/or sign of the QTL effect by the environment. In this context, Multi-Environment Trials (MET), consisting in phenotyping the same panel of genotypes in different locations and/or over different years, have been widely adopted to describe QTL responses to different
environmental factors [25,26]. To detect QTLs of interest from METs, a first approach is to perform single-environment GWAS to identify and further compare QTLs between specific environments [27,28]. A complementary approach is based on multi-environment GWAS in order to detect QTLs with both main effects and QTL × Environments effects (i.e. QTLs with differential effects depending on environmental conditions) [27]. Models have been developed to handle QTL × Environments interactions, such as the multi-trait mixed model MTMM [29] which has been successfully used, firstly in *Arabidopsis thaliana* for the study of flowering time QTLs under different temperatures, and later for plant architecture in barley and in production traits in common bean under contrasted abiotic stress conditions [30–32]. In addition, single-environment GWAS can be exploited via meta-analysis GWAS to identify QTLs with minor effects which might go undetected in separate analyses, as well as to harness markers with stable effects across environments [33,34]. Meta-GWAS are now emerging in plant studies as a practical alternative for combining multiple independent studies to identify additional marker-trait associations and to verify the association previously identified [33,35,36]. A useful aspect of Meta-GWAS for plant is that it can be applied on unbalanced phenotypic datasets (i.e. asymmetrical proportions between class of observation) and it does not require replication of the same individuals across all environments [37].

Taking advantage of the recent publication of reference genomes [38–41] and of the availability of molecular technology such as SNP array [42], GWAS has gradually become more widespread in *Prunus* species for determining the genetic basis of variation in complex agronomic or fruit quality traits [43–47]. Despite these methodological advances, to date G×E interaction effects and the genetic determinants of partial resistance to different biotic stresses from the wide diversity in *Prunus* fruit trees has been little investigated [48–53]. One GWAS experiment has explored non-additive effect in almond and allow the identification of 13 QTLs for several agronomical traits, with only one having additive effects [54]. These results underline that in heterozygous tree crops a significant share of non-additive variance can be exploited in GWAS. Statistical estimates of non-additive effects can therefore allow for a more accurate prediction of total genetic variance which results in a more complete characterization of the genetic architecture of disease resistance traits and help to find complementary QTLs of interest [55,56].

In this context, our objective is to dissect the genetic architecture of pests and diseases resistance or tolerance in peach (*Prunus persica* L.) and apricot (*Prunus armeniaca* L.) in multiple environments by using GWAS. GWAS analysis are particularly well suited for studying peach and apricot genomes by i) being diploid with small size genomes (with respectively 230 and 295 Mbp) [57], ii) harbouring a large untapped phenotypic diversity available in germplasm collection [58,59], and iii) being highly syntenic, thus allowing comparison between them [60]. We built a comprehensive database gathering the
phenotypic response to biotic stresses of peach and apricot using core collections of respectively 206 and 150 accessions monitored under low phytosanitary protection in three and two different orchards in the South-East of France. The pests and diseases observed included leaf curl, leafhopper, powdery mildew, rust and shot hole for peach as well as blossom blight and rust for apricot. We conducted a series of complementary GWAS to i) firstly investigate the main additive and non-additive genetic effects across environments, secondly ii) compare environment specific effects with single-environment GWAS, and finally iii) harness GxE interactions effects using multi-environment GWAS and meta-analysis GWAS.

# Results

## Phenotypic response of peach and apricot accessions to natural infection in multiple orchards

The first step of our analyses aimed at evaluating the biotic performance of the two core collections and to assess the effects of environment, genotype and genotype by environment interactions on biotic stress responses.

### A large variability in the phenotypic response to biotic stress in the two core collections

We scored trees at different time points for each biotic stress throughout the season to fully capture the dynamic expression of the phenotypic diversity under low phytosanitary protection. This allowed us to produce two types of variables for each tree and each biotic stress, namely the i) maximum of susceptibility (Max) and ii) the Area under the Disease Progress Curve (AUDPC).

After spatial correction, we cleaned our dataset in order to retain only environments with sufficient expression of phenotypic variability for each trait, leading to the removal 19.6% of the data (Supplementary Method 2 and Supplementary Fig 1). In the remaining dataset the average standard deviation was of 31.1 and 295.4 respectively for Max and AUDPC (Figure 1 and Supplementary Fig 2). We found the highest diversity of response for apricot rust leaf fall for both Max and AUDPC (with standard deviation of 34.1 and 321.9 respectively). For an equivalent value of Max, trees presented markedly different temporal trajectories for this disease, as illustrated by diverging onset dates and spread of symptoms, with marked environment-specific effects (Figure 1 and Supplementary Fig 3). While correlations between Max and AUDPC remained high for all biotic stresses (from 0.676 to 0.997) (Figure 1 and Supplementary Fig 2), AUDPC thus provided complementary information on the response of accessions to biotic stresses.

To characterize the stability of symptom expression across environment, we calculated the correlations between environment-specific BLUPs for each pest or disease. The genetic correlations between environments ranged from moderate to strong, with average Pearson's correlation
coefficients of 0.68 and 0.73 for Max and AUDPC respectively (Supplementary Table 1 and Supplementary Fig 4). We found the lowest genetic correlation between environments for blossom blight (r = 0.41 in average for Max) and the highest correlations for leaf curl (r = 0.87 and r = 0.91 in average for Max and AUDPC respectively) (Supplementary Table 1). Overall, symptoms were more correlated within a given location than within a given year for both Max and AUDPC (Supplementary Table 1).

## Biotic stress responses are modulated by G and G×E effects

We computed a linear mixed model for each pest and disease to quantify the variance associated to genotype, genotype by location, genotype by year and genotype by year by location interactions (Figure 2.A and Supplementary Method 3). The variance attributed to the genotype amounted between 5.1 and 66.7 % of the total variance. For powdery mildew and shot hole, the genotypic variance was lower than the variance attributed to the G×L interactions (for both Max and AUDPC variables), and for blossom blight lower than the variance attributed to the G×Y interactions (39.7% of total variance). The contribution on the G×L effect on the trait variance was moderate to high (6.3–26.5%) and was higher than the variance explained by G×Y for all biotic stresses, as already expected from previous correlation analysis (Figure 2.A). Finally, G×Y×L variance could not be assessed for all traits and remained moderate (2.7–16.4%). Generally speaking, all responses to biotic stresses were found to be influenced by G×E interactions, with the response to leaf curl being the most stable across environments.

Overall, broad-sense heritability estimates across environments on spatially corrected Max and AUDPC variables ranged from 0.07 to 0.78 for the eight pests and diseases studied (Figure 2.B). Heritability was generally higher when considering AUDPC rather than Max.

When dissecting between additive and non-additive variance, we found that the dominance variance represented up to 16.4% of the total variance, depending on the trait, the highest proportions being found for apricot rust (Supplementary Fig 5). Except for powdery mildew, the dominance variance was higher when considering Max as compared to AUDPC. A comparison of models for these across-environment variance decompositions using the BIC criterion suggested that the Additive-Dominant model was the best model choices for five traits (leafhopper Max and AUDPC, powdery mildew Max and AUDPC, and apricot rust Max).

## Genome-wide association analyses for biotic stress responses in the two core collections

We exploited genetic and G×E variance by applying a series of complementary GWAS models whereby phenotypes were considered i) across environments (G-BLUPs), ii) within environment (environment-
specific BLUPs) and iii) in multiple environments. For each analysis significant, SNPs were grouped in QTLs by considering linkage blocks using the chromosome-specific LD decay values (ranging between 508 and 1620 Kbp and between 131 and 473 bp respectively in the peach and the apricot core collection, depending on the chromosome) (Supplementary Fig 6).

## Complementarity between additive and non-additive QTLs identified across environments

To exploit both additive and non-additive main genetic effects, four different models were tested on G-BLUPs, representing the different modes of action of genes: additive, dominance, recessive and overdominance (the last three have been grouped as 'non-additive'; see material and methods section).

When considering the union of QTLs obtained for Max and AUDPC (i.e. detected for one, the other, or both variables), we identified between 10 and 131 QTLs per trait (Figure 3.A). Within the total of 341 QTLs detected, 188 were additive, 153 were non-additive, and 30 were detected with both additive and non-additive models (Figure 3.A). Except for apricot rust and rust leaf fall, more QTLs were identified with non-additive as with additive models (between 1.2 to 3.7 times more QTLs depending on the trait). Among the non-additive models, the recessive model identified the most QTLs across all traits (between 1 to 28 with an average of 8.8) (Figure 3.C and Supplementary Fig 7), while the over-dominance model identified the least (between 0 to 19 across traits, with an average of 4.9) (Figure 3.D and Supplementary Fig 7). However, it should be reminded that dominant and recessive models are only differentiated regarding the effects of the reference allele. Therefore, QTLs detected with both types of models can correspond to a recessive or dominant mode of action.

We found respectively 271 and 119 QTLs when considering Max and AUDPC variables with a total of 41 overlapping QTLs. Apricot rust was the disease for which Max and AUDPC QTLs overlapped the least (with less than 1% of overlapping QTLs) whereas leaf curl was the disease with the most overlap (with 26% of overlapping QTLs) (Figure 3.A).

Overall, the non-additive GWAS analysis allowed identifying complementary QTLs that were not detected by the additive-GWAS analysis, with 96 and 45 QTLs being exclusive to non-additive models respectively for Max and AUDPC.

## Identification of environment specific QTLs with single-environment GWAS

In order to study environment-specific effects, we implemented single-environment GWAS using single-Locus Mixed-Models on environment-specific BLUPs. We identified between 33 and 568 significant SNPs depending on the trait, corresponding to a total of 801 QTLs (union of results on Max and AUDPC variables) (Figure 4).
The identity and number of associations were greatly dependent on the target environment (62.5% to 100% environment-specific QTLs depending on trait) (Figure 4). Leaf curl was the only trait for which a QTL was identified in all investigated environments and for both AUDPC and Max variables. Apart from this QTL, located on chromosome 4, the other 25 QTLs associated to leaf curl were detected in only one environment. For apricot rust and rust leaf fall, some QTLs have also been identified in all the environments studied when considering AUDPC (representing respectively 5 and 3% of all the QTLs identified). With a total of 6 QTLs detected in 2 different environments and 3 QTLs detected in 3 different environments (total of 5 environments considered), peach rust was the trait with the highest proportion of stable QTLs. On the opposite, blossom blight was the disease with less stable QTLs as all of them were found to be environment-specific. For all traits, QTLs associated to Max appeared to be more stable than those found with AUDPC (22% and 0.09% of QTLs detected in more than one environment, respectively). Finally, among the QTLs identified, between 2 and 46 QTLs overlapped between Max and AUDPC depending on the biotic stress.

Modelling QTLs with GxE interaction effects using MTMM and Meta-GWAS

To statistically distinguish between stable QTLs across environments and those underlying GxE
interactions, we performed i) a multi-environment analysis on environment-specific BLUPs using the
fully parameterized multi-trait mixed model (MTMM)[29] comprising three sub-models ('common
model', 'gei model' and 'full model') (Supplementary Method 6), and ii) a Meta-GWAS with the MetaGE
model [34](Supplementary Method 7). MTMM disentangle common effect markers from interactive
markers and MetaGE jointly analyses single-environment GWAS results to identify QTLs with minor
effects which might go undetected in single environment analyses, as well as markers with significant
effects across multiple environments.

MTMM GWAS and Meta-GWAS analyses led to the identification of respectively 209 and 226 QTLs
significantly associated when considering all traits (union of Max. and AUDPC), with 1 to 61 QTLs per
trait identified with MTMM and 1 to 54 QTLs per trait identified with Meta-GWAS (Figure 5). From
those, 107 QTLs overlapped between MTMM and MetaGE approaches.

MTMM analysis could not be performed on 4 variables (AUDPC of peach rust and leafhopper, and on
both Max and AUDPC for powdery mildew). For the remaining, we identified 75 and 67 QTLs
respectively for ‘common’ and ‘gei’ models (union of Max. and AUDPC, Figure 5). For most traits, more
QTLs have been identified with the ‘common’ model than with the two other MTMM models (between
50 to 100% of QTLs detected). However, for blossom blight (Max), apricot rust (AUDPC and Max) and
apricot rust leaf fall (Max) 53 to 65% of QTLs were identified with the ‘gei’ model. For these ‘gei’,
interactive QTLs, allelic effects of the top SNP ranged from negative to highly positive depending on environment.

When performing Meta-GWAS, we found a large amplitude of z-scores - a summary statistics of the marker's *P*-value of the sign of the marker's effects - for most SNPs. We observed only rare sign inversions, indicative of SNPs with the same qualitative effects (i.e. same direction), but quantitatively different effects across environments (Supplementary Fig 8). A remarkable exception was found for blossom blight, with 61% of the significant SNPs presenting a z-score sign inversion between the 4 environments considered (Supplementary Fig 8.F). Interestingly, we detected 29 QTLs with Meta-GWAS that had not been detected with single environment GWAS (Figure 4 and Figure 5).

## Identification of high confidence QTLs and quantification of their contribution to G and G×E interactions

When considering all traits and models tested, 101 QTLs were detected for peach and 606 for apricot. To retain the most relevant QTLs, high confidence QTLs were defined as those detected in at least 50% of the additive GWAS models tested. We retained up to 32 high confidence QTLs depending on the trait (Supplementary Table 2). Blossom blight was the only disease for which no confidence QTLs could be identified. High confidence QTLs overlapping between AUDPC and Max have been found for leaf curl (1 QTL located on chromosome 4), apricot rust (1 QTL located on chromosome 4) and apricot rust leaf fall (18 QTLs located on chromosome 2, 4 and 6) (Supplementary Table 2).

For each trait, we calculated the proportions of G, G×L and G×Y variance captured by the high confidence QTLs (considered jointly when several confident QTLs were identified for a given pest or disease). The high confidence QTLs captured between 9.6 and 64.6% of the additive genetic variance depending on the biotic stress (Supplementary Fig 9). These QTLs also explained variable portions of G×L interactions (up to 54.4%) and to a lesser extent, of G×Y interactions (up to 29.3%) (Supplementary Fig 9). Interestingly, in apricot rust high confidence QTLs explained a higher proportion of the G×L variance than of the genetic variance (Max variable) (Supplementary Fig 9).

## Local LD and annotation analyses within high confidence QTLs towards identifying candidate genes

Based on the association results, we searched within high confidence QTLs for candidate genes with predicted functions linked to host defence or response mechanisms. As a preliminary step, we readjusted peach QTL intervals to account for the strong LD observed in peach by considering the local LD between the top SNP of each high confidence QTL and all SNPs belonging to this QTL. We considered that a SNP was independent from the top SNP when is r² to top SNP was below 0.2 (Figure 6.A and B
for Leaf_curl_4 and Supplementary Fig 10 for the other peach QTLs). This enabled us to narrow down
the list of candidate genes for peach, with an average of 34.6% of genes with a relevant function
located in the readjusted QTL intervals (Table 1 and Supplementary Table 3).

A total of 87 candidate genes have been identified across the 60 high confidence QTLs, with up to 9 candidate genes for the peach rust QTL on chromosome 2 (Table 1). Among the candidate genes, seven were linked to TMV resistance protein, seven to Leucine-rich receptor-like protein kinase family protein, four to cysteine-rich receptor-like protein kinase, four to TIR-NBS-LRR disease resistance protein and two to MLO-like protein (Table 1). For some QTLs associated to apricot leaf rust, leafhopper and shot hole no candidate genes could be clearly identified.

Interestingly, three apricot candidate genes, 'PruarM.6G017900', 'PruarM.4G204800' and 'PruarM.4G206100', were identified as covering several QTLs (Table 1). We calculated the local LD in these three zones (Figure 6.C. and D, Supplementary Fig 11) and observed two different patterns: for 'PruarM.6G017900' (Figure 6.C) and 'PruarM.4G204800' (Supplementary Fig 11.A) a relatively high LD block overlapped precisely with the underlying candidate gene. These results suggested that the QTLs identified in these different genes were not independent. In contrast, 'PruarM.4G206100' appeared to be located in a moderate LD region with no clear LD blocks (Supplementary Fig 11.B). Max G-BLUPs distributions according to the haplotype constituted by the top SNPs identified in these three genes were plotted (Figure 6.D, for top SNPs belonging to 'PruarM.6G017900' and Supplementary Fig 11.C and D for 'PruarM.4G204800' and "PruarM.4G206100"). For both genes 'PruarM.6G017900' and 'PruarM.4G204800', the two most frequent haplotypes, had significantly different G-BLUPs values (Wilcoxon test P-value < 0.05) (Figure 6.D and Supplementary Fig 11.B).
<table><caption>**Table 1.** List of candidate genes for the high confidence QTLs</caption><thead><tr><th>Biotic stress</th><th>Trait</th><th>QTL Name</th><th>Candidate genes</th><th>Function</th></tr></thead><tbody><tr><td rowspan="5">Leaf curl</td><td rowspan="5">AUDPC and Max</td><td rowspan="5">Leaf_curl_4</td><td>Prupe.4G026800</td><td>LEA hydroxyproline-rich glycoprotein family</td></tr><tr><td>Prupe.4G026900</td><td></td></tr><tr><td>Prupe.4G027000</td><td></td></tr><tr><td>Prupe.4G027300</td><td>Cysteine-rich receptor-like protein kinase</td></tr><tr><td>Prupe.4G027400</td><td></td></tr><tr><td rowspan="3">Leafhopper</td><td rowspan="3">AUDPC</td><td rowspan="3">Leafhopper_7_1</td><td>Prupe.7G065300</td><td>TMV resistance protein N</td></tr><tr><td>Prupe.7G065400</td><td>Disease resistance protein (TIR-NBS-LRR class)</td></tr><tr><td>Prupe.7G065500</td><td></td></tr><tr><td>Powdery mildew</td><td>Max</td><td>Powdery_mildew_6</td><td>Prupe.6G132100</td><td>Disease resistance protein (TIR-NBS-LRR class)</td></tr><tr><td rowspan="10">Peach Rust</td><td rowspan="5">AUDPC</td><td rowspan="5">Peach_Rust_2</td><td>Prupe.2G057600</td><td>TMV resistance protein N</td></tr><tr><td>Prupe.2G057800</td><td></td></tr><tr><td>Prupe.2G058400</td><td>Ethylene-activated signalling pathway</td></tr><tr><td>Prupe.2G058400</td><td></td></tr><tr><td>Prupe.2G059000</td><td></td></tr><tr><td rowspan="3">AUDPC</td><td rowspan="3">Peach_Rust_4</td><td>Prupe.2G059200</td><td>TMV resistance protein N</td></tr><tr><td>Prupe.2G059400</td><td></td></tr><tr><td>Prupe.2G060200</td><td>Ethylene-responsive transcription factor</td></tr><tr><td rowspan="4">AUDPC</td><td rowspan="4">Peach_Rust_4</td><td>Prupe.2G060400</td><td>Disease resistance protein (TIR-NBS-LRR class)</td></tr><tr><td>Prupe.4G051500</td><td>Plant-type hypersensitive response</td></tr><tr><td>Prupe.4G053700</td><td>MLO-like protein</td></tr><tr><td>Prupe.4G054500</td><td>Transcription factor VOZ1</td></tr><tr><td rowspan="2">Shot hole</td><td rowspan="2">AUDPC</td><td rowspan="2">Shot_hole_8</td><td>Prupe.8G257400</td><td>MLO-like protein</td></tr><tr><td>Prupe.8G259700</td><td>Co-chaperone protein</td></tr><tr><td rowspan="10">Apricot Rust</td><td>AUDPC</td><td>Apricot_Rust_2</td><td>PruarM.2G180700</td><td>Leucine-rich receptor-like protein</td></tr><tr><td>AUDPC</td><td>Apricot_Rust_3.2</td><td>PruarM.3G108000</td><td>F-box family protein</td></tr><tr><td>AUDPC</td><td>Apricot_Rust_4.2</td><td>PruarM.4G199900</td><td></td></tr><tr><td>Max</td><td>Apricot_Rust_4.4</td><td></td><td></td></tr><tr><td>AUDPC</td><td>Apricot_Rust_4.5</td><td>PruarM.4G204800</td><td>Leucine-rich receptor-like protein</td></tr><tr><td>AUDPC</td><td>Apricot_Rust_4.6</td><td></td><td></td></tr><tr><td>AUDPC</td><td>Apricot_Rust_4.7</td><td></td><td></td></tr><tr><td>Max</td><td>Apricot_Rust_6.1</td><td>PruarM.6G017900</td><td>Receptor-like serine/threonine-protein kinase RLK1</td></tr><tr><td>Max</td><td>Apricot_Rust_6.2</td><td></td><td></td></tr><tr><td rowspan="7">Apricot Rust Leaf fall</td><td>AUDPC and Max</td><td>Apricot_Rust_leaf_fall_2.1</td><td>PruarM.2G251100</td><td>Disease resistance protein RGA2</td></tr><tr><td>Max</td><td>Apricot_Rust_leaf_fall_2.2</td><td>PruarM.2G283800</td><td>TMV resistance protein N</td></tr><tr><td>Max</td><td>Apricot_Rust_leaf_fall_2.3</td><td>PruarM.2G352900</td><td>WD-40 repeat family protein</td></tr><tr><td>Max</td><td>Apricot_Rust_leaf_fall_2.4</td><td>PruarM.2G359800</td><td>Cysteine-rich receptor-like kinase</td></tr><tr><td>Max</td><td>Apricot_Rust_leaf_fall_2.5</td><td>PruarM.2G359900</td><td></td></tr><tr><td>AUDPC</td><td>Apricot_Rust_leaf_fall_4.1</td><td>PruarM.4G156300</td><td>Lectin protein kinase family protein</td></tr><tr><td>Max and AUDPC</td><td>Apricot_Rust_leaf_fall_4.6</td><td>PruarM.4G199400</td><td></td></tr><tr><td>Max and AUDPC</td><td>Apricot_Rust_leaf_fall_4.7</td><td>PruarM.4G199700</td><td>Leucine-rich receptor-like protein</td></tr></tbody></table>
<table><tr><td>Max and AUDPC</td><td>Apricot_Rust_leaf_fall_4.8</td><td>PruarM.4G199900</td><td></td></tr><tr><td>Max and AUDPC</td><td>Apricot_Rust_leaf_fall_4.10</td><td>PruarM.4G201600</td><td>BURP domain-containing protein</td></tr><tr><td>Max and AUDPC</td><td>Apricot_Rust_leaf_fall_4.16</td><td></td><td></td></tr><tr><td>Max and AUDPC</td><td>Apricot_Rust_leaf_fall_4.17</td><td>PruarM.4G204800</td><td></td></tr><tr><td>Max and AUDPC</td><td>Apricot_Rust_leaf_fall_4.18</td><td></td><td>Leucine-rich receptor-like protein</td></tr><tr><td>Max</td><td>Apricot_Rust_leaf_fall_4.21</td><td>PruarM.4G206100</td><td></td></tr><tr><td>Max and AUDPC</td><td>Apricot_Rust_leaf_fall_4.22</td><td></td><td></td></tr><tr><td>Max</td><td>Apricot_Rust_leaf_fall_6.1</td><td>PruarM.6G017900</td><td>Receptor-like serine/threonine-protein kinase RLK1</td></tr><tr><td>Max and AUDPC</td><td>Apricot_Rust_leaf_fall_6.2</td><td></td><td></td></tr><tr><td>Max</td><td>Apricot_Rust_leaf_fall_6.3</td><td></td><td></td></tr><tr><td>Max and AUDPC</td><td>Apricot_Rust_leaf_fall_6.5</td><td>PruarM.6G018400</td><td>Receptor-like serine/threonine-protein kinase RLK1</td></tr></table>

## Discussion

The main objective of this study was to dissect the genetic architecture of a series of traits reflecting the response to natural infections in multiple environments in peach and apricot. We specifically aimed at deciphering GxE interactions, identifying candidate genomic regions, and comparing the results obtained in these two species.

### Biotic stress response is highly influenced by GxE interactions

Our results revealed that the expression of pest and disease susceptibility was under a strong influence of the environment *per se* but also of GxE interactions (Figure 1, Figure 2, Supplementary Fig 1, Supplementary Fig 2 and Supplementary Fig 3). Variance decompositions (Figure 2.A) and the pairwise genetic correlations (Supplementary Table 1) indicated the prevalence of the GxE interactions over G effects for several traits. Similar patterns have been found in stone fruit diseases, such as in peach brown rot, where skin tolerance to the fungal infection appears to be under high seasonal influence [50,52]. In a study on apricot bacterial canker, Omrani et al. (2019) pointed to the relevance of controlling GxE interactions for accurately estimating genetic effects in the face of environmental variations in multi-year data [49]. In our study, GxE interactions were generally driven more by location than by year effects, which was supported by phenotypic and GWAS analyses. With single-environment GWAS, we identified a total of 567 environment-specific QTLs and 135 environment-shared QTLs with a higher overlap between QTLs identified from each location than from QTLs identified from a given year. Interestingly, the QTLs identified with the multi-environment models
(MTMM and Meta-GWAS) were mostly interactive which changes in direction ('antagonist') or intensity ('differential') according to the environment. Our results thus indicate that data on specific climatic conditions have the potential to reveal additional regions linked to biotic stress response. We thus recommend evaluating the same panel - or highly overlapping panels including several common genotypes - in multiple geographic areas with contrasted environmental conditions.

It should be reminded that the functional validations of G×E interaction QTLs are still rare and requires challenging tests such as gene cloning [61]. Finally, a better characterization of the environments through envirotyping could allow to identify which environmental variables modulate the QTLs effects and to group locations accordingly for common adapted varieties [27].

## A focus on peach leaf curl and apricot blossom blight

Among the seven biotic stresses studied, leaf curl in peach and blossom blight in apricot stood out from the others by presenting two distinct typologies: the former being little influenced by environments and G×E interactions and the latter, by contrast, being strongly influenced by environmental effects. Our study provides the first GWAS analysis on these two major stone fruit diseases.

Apricot blossom blight variation was found to be controlled by a high number of environment-dependent QTLs, indicating a complex genetic determinism. Variance decompositions and GWAS analyses pointed towards predominant G×Y effects on this disease. In fact, all QTLs identified with single-environment models were year-specific and respectively 65% and 61% of the QTLs identified with MTMM and meta-GWAS presented allelic effects modulated by year (i.e. interactive QTLs). Due to these important fluctuations between years, no high-confidence QTLs could be identified with our method. It has been shown that the development of Monilinia on apricot is tightly linked to the interaction between phenology and meteorological conditions, where temperature and humidity during blossom play a major role [62,63], and knowing that the petals being the entry point for the fungus [64,65]. Our dataset was filtered prior to analyses to avoid strong pheno-climatic bias by taking out blossom blight incidence scores for trees having insufficient flower density and Monilinia climatic index in a given year [66]. Despite decent heritability values obtained for apricot blossom blight (0.51; Figure 2.B), we might need to better account for the epidemiological context to further study the genetic basis of that disease.

In contrast to apricot blossom blight, leaf curl damage was a highly heritable trait underlined by low environmental and G×E effects (Figure 2 and Supplementary Table 1). This low environmental impact could be partly explained by a homogeneous leaf curl pressure across the targeted geographic areas, but also within fields and along the susceptibility period (leaf emergence stage). On the opposite of
blossom blight, where the susceptible stage (i.e. flowering) lasts for a few days to a week, the infection period for leaf curl lasts over several weeks [67] meaning that phenology-driven avoidance is rare. In line with the phenotypic analyses, we found a highly stable QTL located on chromosome 4 (identified with all GWAS models tested apart from MTMM). 23 further QTLs have been identified for leaf curl, with an intermediate level of confidence, and mostly driven by location-specific effects (19 location-specific QTLs identified with single-environment models). These results suggest that low susceptibility to leaf curl is governed by a major QTL combined with several minor genes having weaker and environment-specific effects.

In terms of experimental design, we would recommend studying blossom blight incidence over several years in a location with overall high incidence (e.g. upper Rhone Valley in France), to capture for each accession a season during which the tree was actually confronted to the risk. Regarding leaf curl, several locations might be necessary only to capture the minor-effect genes that act as a complement to the major QTL. In terms of breeding strategies, MAS seems to be particularly well adapted to leaf curl, whereas for blossom blight, combining improved epidemiological risk prediction with genomic selection would be a promising strategy.

The rest of the biotic stresses appeared to be more nuanced with no strong typology emerging. We highlighted a location effect but the relative classification of accessions is relatively well preserved between sites (i.e. rank inversions between accessions are not so common). On average, 23% of the QTLs identified with single-environment GWAS were 'environment-shared'. Finally, the majority of 'interactive' QTLs identified presented a 'differential' profile which suggest that some QTLs have effects in some environments and not in others, but true inversions of QTL effects are less common.

**Usefulness of AUDPC and Max variables to study biotic stress response: identification of complementary regions**

We decided to compare Max and AUDPC, two key metrics to capture the complexity of disease dynamics in naturally infected fruit tree orchards. The maximum of damage score provides a clear indicator of 'absolute' vulnerability to biotic stresses but overlooks disease progression over time in different epidemiological contexts.

Although we found high correlations between Max and AUDPC for all biotic stresses studied, with broadly similar heritability (Figure 1, Figure 2.B and Supplementary Fig 2), we identified between 31 to 97 new QTLs with AUDPC depending on the GWAS models. Importantly, among the high confidence QTLs, 32% where only detected using AUDPC. Overall, QTLs detected with Max were more stable while those identified with AUDPC were more location- or year-specific. This underlines the complementary nature of the two metrics. Although widely used in other crops, mostly under controlled infection
conditions, the use of AUDPC has increased in recent years for assessing quantitative disease resistance in fruit trees and for evaluating disease management practices [68–70]. Some studies have even highlighted the dynamics of genetic effects during pathogen colonization with a distinction between early and late response and temporal pattern in effects of QTLs [69].

Because AUDPC requires a large effort to provide repeated measures at multiple time points of the season, we recommend implementing this metric for diseases for which the onset and duration of infection is critical for production, such as apricot rust, or for polycyclic disease [71]. Indeed, early summer defoliation caused by rust has a strong impact on flowering capacity next year. We showed that among the high confidence QTLs identified for apricot rust and rust leaf fall, 25% were AUPDC specific and only 37% overlapped between AUDPC and Max. Therefore, AUDPC enables to better distinguish accessions with a high and continuous level of infection at a critical time point of the season, which is not possible with Max.

## High confidence QTLs reveal the role of basal and host-specific resistance shaping the response to biotic stress in stone fruit trees

Among the QTLs detected with single environment models, respectively 79%, 77% and 54% were also detected with GWAS model focusing on the main genetic effect, by the multi-environment GWAS (MTMM) and with meta-GWAS. In total, these different GWAS models enabled the detection of respectively 38, 79 and 29 QTLs that had not been detected with single environment GWAS. The findings in our study support combining single- and multi-environment models to improve the statistical power, the reliability and robustness of GWAS. Similar results have been obtained in the case of QTL detection for drought response in maize where the incorporation of the interaction with the environment improved the power of GWAS and enabled to find QTL that are significant in a broad range of environments [72]. By conducting GWAS with multiple methods, we narrowed down the list of relevant regions to a total of 60 ‘high confidence QTLs’ across traits. These QTLs explaining in average 36.8% of the genetic variance, and for apricot rust they also explained 53.6% of the genotype by location interaction variance.

To our knowledge, for most of the pests or diseases studied here, resistance or tolerance mechanisms have never been investigated by GWAS analysis neither by QTL mapping, as is the case for leafhopper, shot hole, blossom blight, peach and apricot rust. Previous analyses on peach powdery mildew have already identified QTLs and candidate genes (Vr2 and Vr3) located on chromosomes 8 and 2 [73–75], but none of them overlap with the high confidence QTL identified here. The high confidence QTL we identified for leaf curl on LG4 was not reported previously. A unique linkage analysis study on leaf curl in an interspecific F1 progeny between a peach cultivar and *P. davidiana*, a peach related species,
reported two consistent QTLs located on linkage groups 3 and 6 [76]. QTLs were also identified in these regions by MTMM and Meta-GWAS but they were not categorized as 'high confidence'. In contrast to a QTL-mapping population where the genetic background is homogenous, major alleles might go undetected due to low frequency in a core collection. This may explain the low overlap between our results and those already reported. The type and number of markers included, the approaches used for detecting QTLs or the difference in phenotyping techniques might also contribute to some discrepancies between studies.

Interestingly, the majority of candidate genes identified belong to the Leucine-rich repeat (LRR)-containing receptors (LRR-CRs) family gene which represents the largest class of plant resistance genes (R genes) [77,78]. LRR-CRs can be classified in three main sub-families: LRR receptor-like kinase (LRR-RLK), LRR receptor-like protein (LRR-RLP) and nucleotide-binding-site LRR (NBS-LRR) [79,80]. In our study we identified candidate genes belonging to all the categories with genes related to i) LRR-RLK identified for apricot rust and rust leaf fall (PruarM.6G017900 and PruarM.6G018400), ii) LRR-RLP identified for apricot rust and rust leaf (PruarM.2G180700, PruarM.4G204800, PruarM.4G199400, PruarM.4G199700, PruarM.4G199900, PruarM.4G204800 and PruarM.4G206100) and iii) NBS-LRR identified, for leafhopper (Prupe.7G065400 and Prupe.7G065500), powdery mildew (Prupe.6G132100), peach rust (Prupe.2G060400). RLKs and RLPs have been identified to function as cell surface pattern-recognition receptors (PRRs) which detect pathogens-associated molecular patterns (PAMPs) and lead to the activation of the pattern-triggered immunity (PTI), the basal defence responses to pathogens [79,81]. We also identified a high number of genes annotated as 'TMV resistance protein N' (Prupe.7G065300 identified for leafhopper, Prupe.2G057600, Prupe.2G057800, Prupe.2G059000, Prupe.2G059200 and Prupe.2G059400 identified for peach rust, and PruarM.2G283800 for apricot rust leaf fall). The tobacco N gene, conferring resistance to tobacco mosaic virus (TMV) has been identified as TIR-NBS-LRR protein type [82]. It has been shown that NBS-LRR proteins conferred resistance to diverse pathogens, including fungi, oomycetes, bacteria, viruses and insects. By recognizing pathogen effector molecules, NBS-LRR proteins lead to the activation of defence responses and play a fundamental role in the effector-triggered immunity (ETI), the second pathogen-sensing mechanism in plants [81]. ETI usually induces programmed cell death at the infection site, known as the hypersensitive response (HR), and thus locally limits pathogen spread. NBS-LRR proteins are highly variable and subject to balancing selection, which could contribute to the high specificity of these effectors [79]. Thus, our results point towards a prevalent role of host specific (narrow spectrum) and non-specific (large spectrum) resistance mechanisms among the diseases studied here.
Similar, yet divergent: difference in the architecture of biotic stress response between peach and apricot in the light of their biological features

Peach and apricot are closely related species belonging to the *Rosaceae* family and the genus *Prunus*. They present different reproductive characteristics and demographic history [83]. Endowed with a self-incompatibly system, apricot is a cross-pollinating species [84], whereas peach is a preferentially self-pollinating species [85]. This results for apricot in a high number of effective recombination events and a rapid expected LD decay over the genome, whereas the opposite pattern is expected for peach [86]. In our study, global mean LD blocs over all chromosomes are about 3300 times longer for the peach core collection compared to the apricot core collection. In fact, previous association studies have already highlighted a LD decay over very short distances (100 to 200 bp) in apricot (based on a population comprising an important part of the material used in presented work) [48,49]. In a similar way for peach, other studies have found a significant LD decay at 5 Kbp to 1.8 Mbp (depending on the population), confirming high levels of LD over long distance for that crop [44,87]. QTL intervals were therefore particularly small for apricot (max 946 Bp), with most of them being defined by the top SNP. Minding a sufficient number of markers, a higher mapping resolution makes association analyses more powerful to target causal polymorphisms in apricot, facilitating the identification of promising candidate genes [88]. For rust tolerance, we identified three genes covering several confidence QTLs in apricot, along with favourable haplotypes (Figure 6.D and Supplementary Fig 11). Results on apricot strongly contrast with those obtained in peach, where QTLs were up to max 3240 Kbp in size, and contained 52 to 528 genes. Using local LD values allowed to refine QTL sizes and to narrow down the list to 23 candidate genes.

In our study, rust was the only disease able to affect both peach and apricot. Rust diseases have a fungal origin and affect many crops, and typically present a high specificity to their host [89]. *Tranzchelia discolor* (Fuckel) has been shown to attack peach, almond, plum, apricot, and cherry [90,91]. For both species, a high confidence QTL was identified on chromosome 2.

Finally, we identified a total of 123 new QTLs using non-additive models, indicating that non-additive effects play a significant role in biotic stress response for both species. In proportion, more non-additive QTLs were found in peach (56%) than in apricot (39%). Although peach is a self-compatible species with a low level of heterozygosity, this result indicates that both additive and non-additive genetic effects are important in the genetic architecture of the traits under study, as already been shown for leaf curl susceptibility inheritance [92]. In almond, another cross-pollinating *Prunus* species, non-additive GWAS conducted on agronomic traits led to the identification of 13 QTLs with only one
having additive effect [83]. Harnessing non-additive variance might thus be valuable to better exploit genetic variance, which is particularly relevant in clonally propagated crops.

## Conclusion and prospect

This study provides a thorough characterization of the genetic variation underlying resistance and tolerance to several pests or diseases in peach and apricot via observations over multiple environments under natural conditions. By quantifying and exploiting G×E interactions, we found a majority of environment-specific or interactive QTLs. For future experiments, we recommend considering the combination of pathogens which are specifically expressed in each target location to optimize the choice of experimental sites.

Our results emphasize that, while accounting for pedo-climatic adaptation and market requirements, future breeding should focus on adaptation to the biotic profile of the target production area. Developing new peach and apricot cultivars is a time-consuming process because of the long juvenile stage, large plant size, and the multi-year recording of phenotypic performance [93]. As a next step, the effects of the allelic variants identified in this study need to be validated before incorporating them in MAS strategy to efficiently screen breeding materials and optimize the multi-trait and multi-environment selection schemes. Another promising strategy to overcome the limiting aspects of GWAS would be the use of genomic prediction models for a better exploitation of weak effects underlying quantitative resistances and tolerances [94,95].

## Material and methods

### Plant material and trial description

The plant material and experimental design used in this study were previously described in [96]. To sum up, we used two core collections of respectively 206 peach (P. persica and 3 relatives: P. davidiana, P. mira, and P. kansuensis) and 150 apricot (P. armeniaca) accessions coming from multiple geographical origins with significant number of accessions of historical and regional importance, landraces and elite materials. The peach core collection also includes ornamental cultivars, few accessions from close-related species and interspecific crosses.

All orchards are located in South-East France, in three locations for peach: Étoile-sur-Rhône (44.8°N/ 4.9°E, Drôme), Avignon (43.9°N/4.8°E, Vaucluse) and Torreilles (42.8°N/3.0°E, Pyrénées-Orientales), respectively called ET, AV and TO; and in two different locations for apricot: Saint-Marcel-lès-Valence (45.0°N/4.9°E, Drôme) and Bellegarde (43.7°N/4.5°E, Gard), respectively called SM and BE. Climate ranges from Mediterranean to semi-continental. Peach and apricot orchards are composed of two and five randomized complete blocks respectively, each block containing one replicate per genotype. All
accessions have been grafted on 'Montclar® Chanturgue' peach rootstocks. First leaf emergence was in 2019 for peach and in 2018 for apricot.

The orchards were managed under integrated plant protection during the first two years after the plantation with a progressive reduction of phytosanitary protection thereafter. Overall, only punctual treatments were applied if the survival of the orchard was at stake for peach orchards since 2023, and except for the protection against the psyllid vector of the Apricot Chlorotic Leaf Roll (ACLR) in Bellegarde, no protection was applied for apricot orchards since 2020.

## Phenotypic data collection

We monitored the damages caused by four diseases and one pest on peach in 2021, 2022 and 2023 and by two diseases on apricot in 2020, 2021, 2022 and 2023 (Table 2). For the rest of the analysis, the term of environment refers to the combination of a location and a year. All field data were collected using the 'Field Book' Android application [97].

The assessment of the damages caused by all these pests and diseases was based on the visual estimation of symptom incidence (percentage of damaged leaves). Different rating scales were used to best capture the phenotypic expression of each pest or disease (Supplementary Method 1).

Except for blossom blight which was assessed around 30 days after blooming date and for which the shoots infected by *Monilinia* spp. were removed, repeated scoring was undertaken for the different pests and diseases to monitor the evolution of contamination throughout the seasons in the different orchards.

## Genotyping and resequencing

Genotyping of the peach core collection was performed on 192 individuals using the higher density IRSC 16K SNP array (https://www.rosaceae.org/Analysis/431). The genotyping dataset was filtered to retain SNPs with call rate per marker > 90%, missingness per individual < 50%, heterozygosity ≤ 98% and MAF > 5% which resulted in a final set 13,680 markers.

The apricot core collection was sequenced with the Illumina HiSeq 2000 NGS technique. The sequence alignment was performed on the third version of the 'Marouch' genome [41]. The dataset was first filtered with a missingness per individual < 50% and only biallelic SNPs with more than 10 reads deep/SNP have been conserved. Finally, the dataset was filtered to retain SNPs with call rate per marker > 95%, heterozygosity ≤ 98%, MAF > 5% and a removal of duplicated SNPs. This dataset consisted in 207,000 markers.

For both datasets, the imputation of missing data was done using the most frequent allele.
For the datasets used in non-additive GWAS, we included two more filtering steps to retain SNPs with three genotypic classes and with a minimum genotypic class frequency higher than 5%. These datasets consisted in respectively 10,082 markers for peach and 25,508 markers for apricot.

<table><caption>**Table 2.** Summary of pests and diseases monitored in this study</caption><thead><tr><th></th><th>Pest and Disease</th><th>Organism</th><th>Pathogen</th><th>Organ</th><th>Number of environments observed</th><th>Number of environments considered</th></tr></thead><tbody><tr><td rowspan="5">Peach</td><td>Leaf curl</td><td>Fungi</td><td>Taphrina deformans</td><td>Leaves</td><td>3 (3 locations and 1 year)</td><td>3 (3 locations and 1 year)</td></tr><tr><td>Leafhopper</td><td>Insect</td><td>Asymmetrasca decedens</td><td>Leaves</td><td>6 (3 locations and 3 years)</td><td>6 (3 locations and 3 years)</td></tr><tr><td>Powdery mildew</td><td>Fungi</td><td>Podosphaera pannosa</td><td>Leaves</td><td>8 (3 locations and 3 years)</td><td>5 (2 locations and 3 years)</td></tr><tr><td>Rust</td><td>Fungi</td><td>Tranzchelia spp.</td><td>Leaves</td><td>6 (2 locations and 3 years)</td><td>5 (2 locations and 3 years)</td></tr><tr><td>Shot hole</td><td>Fungi</td><td>Wilsonomyces carpophilus</td><td>Leaves</td><td>9 (3 locations and 3 years)</td><td>4 (3 locations and 2 years)</td></tr><tr><td rowspan="2">Apricot</td><td>Blossom and twig blight</td><td>Fungi</td><td>Monilinia spp. (mostly Monilia laxa)</td><td>Flowers and twigs</td><td>4 (1 location and 4 years)</td><td>4 (1 location and 4 years)</td></tr><tr><td>Rust</td><td>Fungi</td><td>Tranzschelia spp.</td><td>Leaves</td><td>6 (2 locations and 4 years)</td><td>6 (2 locations and 4 years)</td></tr></tbody></table>

## Statistical analysis of phenotypic data

All statistical analyses of phenotypic and genotypic data were performed using the R software version 4.1.2 (R-Core-Team, 2020)

At each date of observation, phenotypic data were adjusted within each environment to correct for spatial environmental effects using information on row and column (Supplementary Method 2). After removing environments with insufficient phenotypic expression, we extracted Max and AUDPC (Supplementary Method 2). These variables were transformed using Box-Cox to correct for heteroscedasticity and non-normality of error terms (Supplementary Method 2). Variance analysis was performed using mixed linear models to quantify the variance associated to genetic, GxE, and non-additive effects (Supplementary Method 3 and Supplementary Method 4).

## Linkage disequilibrium

Intrachromosomic linkage disequilibrium (LD) was measured with the r² value estimated by using VCFTools v0.1.16 [98]. The r² was calculated for every combination of SNPs within a sliding window of
10 Mbp for peach and of 30 Kbp for apricot for every chromosome individually. We used a threshold of 0.2 to set the LD decay which was then represented graphically using Hill and Weir equation [99].

Genome-wide association studies (GWAS)

<table><caption>**Table 3.** Summary of the different GWAS model used in this study</caption><thead><tr><td>GWAS name</td><td>R Package</td><td>Rationale</td></tr></thead><tbody><tr><td>Main genetic effect additive GWAS</td><td>MM4LMM</td><td>Decipher the main additive genetic effects across environments</td></tr><tr><td>Main genetic effect non-additive GWAS</td><td>MM4LMM</td><td>Decipher the main non-additive genetic effect across environments</td></tr><tr><td>Single environment GWAS</td><td>MM4LMM</td><td>Study genetic effects environment by environment</td></tr><tr><td>Multi-environment GWAS</td><td>MTMM</td><td>Assess the effect of G × E interactions on QTL detection and identify genetic factors with common effects, i.e. stable marker effects across environments, or interaction effects, i.e. differential marker effects across environments</td></tr><tr><td>Meta-analysis GWAS</td><td>MetaGE</td><td>Assess the effect of G × E interactions on QTL detection by a joint analysis of single-environment GWAS results</td></tr></tbody></table>

Across-environment GWAS based on main additive and non-additive genetic as well as single-environment GWAS (Table 3) were performed using the single locus linear mixed model of [100] implemented in the R package MM4LMM [101]:

$$Y = X\beta + Qv + Zu + \varepsilon \quad (1)$$

where *Y* is respectively the vector of G-BLUPs for across-environment GWAS or of environment-specific BLUPs for single-environment GWAS, *X* is the molecular marker score matrix, *β* is the vector of marker effects, *Q* is the structure matrix estimated by the SNMF R package [102] (Supplementary Method 5 and Supplementary Result 1), *v* is the fixed effect of the structure, *Z* is an incidence matrix, *u* is the vector of random background polygenic effects with variance $\sigma^2_u = K \sigma^2_G$ (where *K* is the variance-covariance matrix of *G*, and $\sigma^2_G$ is the genetic variance), and $\varepsilon$ is the vector of residuals.

The estimation of the variance-covariance matrix of *G* was determined by a genetic relatedness (or kinship) matrix, derived from all SNPs except those on the chromosome containing the SNP being tested as in [103]. According to the model, additive or non-additive kinships were used. They were calculated using the method developed by [104] and the source code published by [105]. For every pest or disease, we tested two different corrections including the kinship (model 'K') or both kinship and population structure (model 'Q + K'). The best model was chosen on a trait-by-trait basis by comparing the likelihood of each null model (i.e. with only kinship of kinship + structure effects) using
the AIC. We used a 'modified' 5% Bonferroni correction described by [106] to define significant associations between phenotypic data and genotypic markers. This threshold is based on the effective number of independent SNPs calculated thanks to the *Meff* equation implemented in the R package *poolR* [107]. The number of effective markers was 1326 for peach and 1778 for apricot.

To explore non-additive genotype-phenotype associations, we transformed the genotypic datasets following [108] and [54], as illustrated in Supplementary Table 4.

Subsequently, we performed single locus multi-environment GWAS by using the multi-trait mixed model (MTMM) as described by [29] in order to distinguish common and environment-specific marker-phenotype association (Table 3 and Supplementary Method 6). We finally performed a meta-analysis relying on the results of the single-environment GWAS model by using MetaGE R package [34] (Table 3 and Supplementary Method 7).

## Identification of high confidence QTLs and variance analysis of GxE and QTLxE

For all GWAS analyses, every QTL region was defined using the position of the top SNP (i.e. with the lowest *P*-value) and the estimated LD decay for every chromosome. In other words, QTLs were centred on the top SNP position with an extent of 2 time-specific chromosomal LD decay. Candidate SNPs distant less than the LD decays size were considered as belonging to this QTL. Each QTL name was defined as the associated trait name, followed by the chromosome number and the QTL's order among all QTLs for this trait along the chromosome.

For every trait, we selected 'high confidence' QTLs which correspond to QTLs significant in at least 50% of the different additive GWAS model tested for the maximum of damage score or for AUDPC. We considered that each environment of the single-environment GWAS as separate models as well as the model 'full', 'common' and 'gei' of MTMM.

To assess the biological relevance of the high confident QTLs, we followed the approach previously described by [27,109] which aim at comparing the variance of a multi-locus multi-environment mixed model with and without the fixed effects of the QTLs and of the QTLs by environment interaction (Supplementary Method 8). Last, we calculated the indicator $\gamma_{q,r}$ to assess the proportion of genetic, G×L interaction and G×Y interaction variance explained jointly by all detected QTLs for each trait (Supplementary Method 8).

## Candidate gene analysis

Each high confidence QTLs region was queried in the *Prunus persica* Whole Genome v2.0 Assembly and Annotation v2.1 [39] or in the *Prunus armeniaca* Marouch whole genome v1.0 [41] from the
Genome Database for Rosaceae (www.rosaceae.org). SNP gene localization and function were determined using the JBrowse tool (www.rosaceae.org/tools/jbrowse). Their function was then investigated in the literature, and genes involved in traits linked to host defence mechanisms were therefore proposed as putative candidate genes.

For peach, the local LD was calculated between the top SNP of the different ‘high confidence’ QTLs and all the SNPs belonging in this QTLs by using LDcorSV R package which takes into account population stratification and the relatedness between genotypes [110]. We used r² > 0.2 to filter the putative candidate genes.

For genes which overlapped among several confidence QTLs for apricot, pairwise LD from all physically close SNPs in the region was computed using r² calculated with LDcorSV R package [110]. The haplotypes were constructed by considering all top SNPs of the confidence QTLs involved for each G-BLUPs phenotype. Mean phenotypic distributions of the different haplotypes were compared using paired samples Wilcoxon tests for the associated G-BLUPs phenotypes.

## Acknowledgements and Funding

This work was supported by funding through LabEx AGRO 2011-LABX-002 (under I-Site Muse framework) coordinated by Agropolis Foundation (project ID: 2002-030), the INRAE Department for Plant Genetics and Breeding, the France AgriMer CASDAR Project ‘RésiDiv’ (project ID: 6846752) and the EU Horizon Innovation Actions InnOBreed n°101061028.

We are thankful to the two experimental research stations SEFRA and SICA-Centrex, to the INRAE experimental unit A2M and to the INRAE research experimental unit UERI of Gotheron for the maintenance of the two core collections. We thank Thierry Pascal for initiating the peach collection. We also thank Freddy Combe (INRAE, UERI Gotheron), Amandine Fleury (INRAE, UERI Gotheron), Luana Gillet (INRAE, UR GAFL), Béatrice Monnet (INRAE, UE A2M) as well as several students for their valuable contribution in field experiment and plant phenotyping. We are also grateful to Guillaume Roch (CEP Innovation) and Guy Clauzel (CEP Innovation) for grafting, pruning and DNA sampling the apricot core collection, to Jacques Lagnel (INRAE, UR GAFL) and Frédérique Bitton (INRAE, UR GAFL) for coordinating the bioinformatic processing of the raw data from the re-sequencing of the apricot core collection and to Océane Perez (INRAE, UR GAFL) for the preliminary analysis on apricot. We acknowledge Emilie Millet (INRAE, UR GAFL) for sharing her expertise in statistical analysis and providing valuable scientific advice.
## Author contributions

MR, BQ and JMA conceived the project and designed the experiment. JMA, BQ, AB and LB created the core collections. ND collected the DNA samples and coordinated the genotyping of the peach core collection. JMA collected the DNA samples and coordinated the resequencing of the apricot core collection. LH contributed to the alignment and filtration of the raw data from re-sequencing. MS, AB, LB, FG, MLP, MR, VSi and SV realised the phenotypic observations. MS analysed the data under the supervision of MR. VSe supervised GWAS analyses. LH contributed to the identification of candidate genes. MS wrote the manuscript. MR, BQ, JMA and VSe reviewed and edited the manuscript. All authors approved the manuscript.

## Data availability statement

The raw phenotypic and genotypic datasets and the R scripts generated for this study are accessible on the data.gouv public repository at: https://doi.org/10.57745/2HNRN0

## Conflict of interest

The authors declare no conflict of interest.
## References

1. Oerke E-C. Crop losses to pests. *The Journal of Agricultural Science* 2006;**144**:31–43.

2. Savary S, Bregaglio S, Willocquet L *et al*. Crop health and its global impacts on the components of food security. *Food Security* 2017;**9**:311–27.

3. Deng Y, Ning Y, Yang D-L *et al*. Molecular Basis of Disease Resistance and Perspectives on Breeding Strategies for Resistance Improvement in Crops. *Molecular Plant* 2020;**13**:1402–19.

4. Agrios GN. *Plant Pathology*. Elsevier, 2005.

5. Råberg L, Graham AL, Read AF. Decomposing health: tolerance and resistance to parasites in animals. *Philosophical Transactions of the Royal Society B: Biological Sciences* 2008;**364**:37–49.

6. Mundt CC. Durable resistance: A key to sustainable management of pathogens and pests. *Infection, Genetics and Evolution* 2014;**27**:446–55.

7. Gessler C, Patocchi A, Sansavini S *et al*. Venturia inaequalis Resistance in Apple. *Critical Reviews in Plant Sciences* 2006;**25**:473–503.

8. Kleinhofs A, Brueggeman R, Nirmala J *et al*. Barley Stem Rust Resistance Genes: Structure and Function. *The Plant Genome* 2009;**2**, DOI: 10.3835/plantgenome2009.02.0011.

9. Palloix A, Ayme V, Moury B. Durability of plant major resistance genes to pathogens depends on the genetic background, experimental evidence and consequences for breeding strategies. *New Phytologist* 2009;**183**:190–9.

10. Brun H, Chèvre A-M, Fitt BD *et al*. Quantitative resistance increases the durability of qualitative resistance to Leptosphaeria maculans in Brassica napus. *New Phytologist* 2010;**185**:285–99.

11. Baumgartner IO, Patocchi A, Frey JE *et al*. Breeding Elite Lines of Apple Carrying Pyramided Homozygous Resistance Genes Against Apple Scab and Resistance Against Powdery Mildew and Fire Blight. *Plant Mol Biol Rep* 2015;**33**:1573–83.

12. Wiesner-Hanks T, Nelson R. Multiple disease resistance in plants. *Annual Review of Phytopathology* 2016;**54**:229–52.

13. Foulongne M, Pascal T, Pfeiffer F *et al*. QTLs for powdery mildew resistance in peach ×Prunus davidiana crosses: consistency across generations and environments. *Molecular Breeding* 2003;**12**:33–50.

14. Sauge M-H, Lambert P, Pascal T. Co-localisation of host plant resistance QTLs affecting the performance and feeding behaviour of the aphid Myzus persicae in the peach tree. *Heredity* 2012;**108**:292–301.

15. Cirilli M, Rossini L, Geuna F *et al*. Genetic dissection of Sharka disease tolerance in peach (P. persica L. Batsch). *BMC Plant Biology* 2017;**17**:1–15.

16. Hajjar R, Hodgkin T. The use of wild relatives in crop improvement: a survey of developments over the last 20 years. *Euphytica* 2007;**156**:1–13.
17. Tibbs Cortes L, Zhang Z, Yu J. Status and prospects of genome-wide association studies in plants. *Plant Genome* 2021; **14**, DOI: 10.1002/tpg2.20077.

18. Desaint H, Hereil A, Causse M. Genome-Wide Association Study: A Powerful Approach to Map QTLs in Crop Plants. In: Raina A, Wani MR, Laskar RA, et al. (eds.). *Advanced Crop Improvement, Volume 1*. Cham: Springer International Publishing, 2023, 379–455.

19. Monnot S, Desaint H, Mary-Huard T *et al*. Deciphering the Genetic Architecture of Plant Virus Resistance by GWAS, State of the Art and Potential Advances. *Cells* 2021; **10**:3080.

20. Xu Y, Crouch JH. Marker-Assisted Selection in Plant Breeding: From Publications to Practice. *Crop Science* 2008; **48**:391–407.

21. Hasan N, Choudhary S, Naaz N *et al*. Recent advancements in molecular marker-assisted selection and applications in plant breeding programmes. *Journal of Genetic Engineering and Biotechnology* 2021; **19**:128.

22. Sibiya J, Tongoona P, Derera J *et al*. Genetic analysis and genotype × environment (G × E) for grey leaf spot disease resistance in elite African maize (Zea mays L.) germplasm. *Euphytica* 2012; **185**:349–62.

23. Lussewa RK, Edema R, Lamo J. Magnitude of genotype x environment interaction for bacterial leaf blight resistance in rice growing areas of Uganda. *African Crop Science Journal* 2016; **24**:11–24.

24. Zeng Y, Shi J, Ji Z *et al*. Genotype by Environment Interaction: The Greatest Obstacle in Precise Determination of Rice Sheath Blight Resistance in the Field. *Plant Disease* 2017; **101**:1795–801.

25. Boer MP, Wright D, Feng L *et al*. A Mixed-Model Quantitative Trait Loci (QTL) Analysis for Multiple-Environment Trial Data Using Environmental Covariables for QTL-by-Environment Interactions, With an Example in Maize. *Genetics* 2007; **177**:1801–13.

26. Rodrigues PC. An overview of statistical methods to detect and understand genotype-by-environment interaction and QTL-by-environment interaction. *Biometrical Letters* 2018; **55**:123–38.

27. Millet EJ, Welcker C, Kruijer W *et al*. Genome-Wide Analysis of Yield in Europe: Allelic Effects Vary with Drought and Heat Scenarios. *Plant Physiology* 2016; **172**:749–64.

28. Hu X, Wang G, Du X *et al*. QTL analysis across multiple environments reveals promising chromosome regions associated with yield-related traits in maize under drought conditions. *The Crop Journal* 2021; **9**:759–66.

29. Korte A, Vilhjálmsson BJ, Segura V *et al*. A mixed-model approach for genome-wide association studies of correlated traits in structured populations. *Nat Genet* 2012; **44**:1066–71.

30. Sasaki E, Zhang P, Atwell S *et al*. “Missing” G × E Variation Controls Flowering Time in Arabidopsis thaliana. *PLOS Genetics* 2015; **11**:e1005597.

31. Oladzad A, Porch T, Rosas JC *et al*. Single and Multi-trait GWAS Identify Genetic Factors Associated with Production Traits in Common Bean Under Abiotic Stress Environments. *G3 Genes / Genomes / Genetics* 2019; **9**:1881–92.

32. Bretani G, Shaaf S, Tondelli A *et al*. Multi-environment genome-wide association mapping of culm morphology traits in barley. *Front Plant Sci* 2022; **13**, DOI: 10.3389/fpls.2022.926277.
33. Shook JM, Zhang J, Jones SE et al. Meta-GWAS for quantitative trait loci identification in soybean. *G3 Genes/Genomes/Genetics* 2021;**11**:jkab117.

34. De Walsche A, Vergne A, Rincent R et al. *metaGE: Investigating Genotype × Environment Interactions through Meta-Analysis*. Plant Biology, 2023.

35. Zhao J, Sauvage C, Zhao J et al. Meta-analysis of genome-wide association studies provides insights into genetic control of tomato flavor. *Nat Commun* 2019;**10**:1534.

36. Joukhadar R, Thistlethwaite R, Trethowan R et al. Meta-analysis of genome-wide association studies reveal common loci controlling agronomic and quality traits in a wide range of normal and heat stressed environments. *Theor Appl Genet* 2021;**134**:2113-27.

37. Evangelou E, Ioannidis JPA. Meta-analysis methods for genome-wide association studies and beyond. *Nat Rev Genet* 2013;**14**:379-89.

38. Verde I, Abbott AG, Scalabrin S et al. The high-quality draft genome of peach (Prunus persica) identifies unique patterns of genetic diversity, domestication and genome evolution. *Nature Genetics* 2013;**45**:487-94.

39. Verde I, Jenkins J, Dondini L et al. The Peach v2.0 release: high-resolution linkage mapping and deep resequencing improve chromosome-scale assembly and contiguity. *BMC Genomics* 2017;**18**:1-18.

40. Jiang F, Zhang J, Wang S et al. The apricot (Prunus armeniaca L.) genome elucidates Rosaceae evolution and beta-carotenoid synthesis. *Horticulture Research* 2019;**6**:128.

41. Groppi A, Liu S, Cornille A et al. Population genomics of apricots unravels domestication history and adaptive events. *Nat Commun* 2021;**12**:3956.

42. Mas-Gómez J, Cantín CM, Moreno MÁ et al. Exploring Genome-Wide Diversity in the National Peach (Prunus persica) Germplasm Collection at CITA (Zaragoza, Spain). *Agronomy* 2021;**11**:481.

43. Dhanapal AP, Crisosto CH. Association genetics of chilling injury susceptibility in peach (Prunus persica (L.) Batsch) across multiple years. *3 Biotech* 2013;**3**:481-90.

44. Micheletti D, Dettori MT, Micali S et al. Whole-Genome Analysis of Diversity and SNP-Major Gene Association in Peach Germplasm. Han Y (ed.). *PLoS ONE* 2015;**10**:e0136803.

45. Cao K, Zhou Z, Wang Q et al. Genome-wide association study of 12 agronomic traits in peach. *Nat Commun* 2016;**7**:13246.

46. Font i Forcada C, Guajardo V, Chin-Wo SR et al. Association Mapping Analysis for Fruit Quality Traits in Prunus persica Using SNP Markers. *Front Plant Sci* 2019;**9**, DOI: 10.3389/fpls.2018.02005.

47. Tan Q, Li S, Zhang Y et al. Chromosome-level genome assemblies of five Prunus species and genome-wide association studies for key agronomic traits in peach. *Hortic Res* 2021;**8**:213.

48. Mariette S, Wong Jun Tai F, Roch G et al. Genome-wide association links candidate genes to resistance to Plum Pox Virus in apricot (Prunus armeniaca). *New Phytologist* 2016;**209**:773-84.
49. Omrani M, Roth M, Roch G et al. Genome-wide association multi-locus and multi-variate linear mixed models reveal two linked loci with major effects on partial resistance of apricot to bacterial canker. *BMC Plant Biology* 2019;**19**:1–18.

50. Fu W, da Silva Linge C, Gasic K. Genome-Wide Association Study of Brown Rot (Monilinia spp.) Tolerance in Peach. *Front Plant Sci* 2021;**12**, DOI: 10.3389/fpls.2021.635914.

51. Li X, Wang J, Su M et al. Single Nucleotide Polymorphism Detection for Peach Gummosis Disease Resistance by Genome-Wide Association Study. *Front Plant Sci* 2022;**12**:763618.

52. Martínez-García PJ, Mas-Gómez J, Prudencio ÁS et al. Genome-wide association analysis of Monilinia fructicola lesion in a collection of Spanish peach landraces. *Frontiers in Plant Science* 2023;**14**, DOI: 10.3389/fpls.2023.1165847.

53. Li X, Wang J, Su M et al. Multiple-statistical genome-wide association analysis and genomic prediction of fruit aroma and agronomic traits in peaches. *Horticulture Research* 2023;**10**:uhad117.

54. Pérez de los Cobos F, Coindre E, Dlalah N et al. Almond population genomics and non-additive GWAS reveal new insights into almond dissemination history and candidate genes for nut traits and blooming time. *Horticulture Research* 2023;**10**:uhad193.

55. Tsepilov YA, Shin S-Y, Soranzo N et al. Nonadditive Effects of Genes in Human Metabolomics. *Genetics* 2015;**200**:707–18.

56. Xue Y, Liu S, Li W et al. Genome-Wide Association Study Reveals Additive and Non-Additive Effects on Growth Traits in Duroc Pigs. *Genes* 2022;**13**:1454.

57. Arumuganathan K, Earle ED. Nuclear DNA content of some important plant species. *Plant Molecular Biology Reporter* 1991;**9**:208–18.

58. Li X, Meng X, Jia H et al. Peach genetic resources: diversity, population structure and linkage disequilibrium. *BMC Genet* 2013;**14**:84.

59. Bourguiba H, Scotti I, Sauvage C et al. Genetic Structure of a Worldwide Germplasm Collection of Prunus armeniaca L. Reveals Three Major Diffusion Routes for Varieties Coming From the Species' Center of Origin. *Frontiers in Plant Science* 2020;**11**, DOI: 10.3389/fpls.2020.00638.

60. Dirlewanger E, Graziano E, Joobeur T et al. Comparative mapping and marker-assisted selection in Rosaceae fruit crops. *Proc Natl Acad Sci USA* 2004;**101**:9891–6.

61. Marais DLD, Hernandez KM, Juenger TE. Genotype-by-Environment Interaction and Plasticity: Exploring Genomic Responses of Plants to the Abiotic Environment. *Annual Review of Ecology, Evolution, and Systematics* 2013;**44**:5–29.

62. Tamm. Phenological Analysis of Brown Rot Blossom Blight of Sweet Cherry Caused by Monilinia laxa. *Phytopathology* 1995, DOI: 10.1094/Phyto-85-401.

63. Luo Y, Morgan DP, Michailides TJ. Risk analysis of brown rot blossom blight of prune caused by Monilinia fructicola. *Phytopathology* 2001;**91**:759–68.

64. Holb IJ. Brown rot blossom blight of pome and stone fruits: symptom, disease cycle, host resistance, and biological control. *International Journal of Horticultural Science* 2008;**14**:15-21.
65. Brun L, Tresson P, Dubois J et al. Suivre la contamination des fleurs par Monilinia laxa—Le développement d'un modèle phéno-climatique de contamination des fleurs d'abricotier par la moniliose va permettre d'optimiser le pilotage de la protection contre la maladie. *Phytoma* 2021;740:38-42.

66. Tresson P, Brun L, García de Cortázar-Atauri I et al. Future development of apricot blossom blight under climate change in Southern France. *European Journal of Agronomy* 2020;112:125960.

67. Giosuè S, Spada G, Rossi V et al. Forecasting Infections of the Leaf Curl Disease on Peaches Caused by Taphrina deformans. *European Journal of Plant Pathology* 2000;106:563-71.

68. Osorio-Guarín JA, Berdugo-Cely JA, Coronado-Silva RA et al. Genome-Wide Association Study Reveals Novel Candidate Genes Associated with Productivity and Disease Resistance to Moniliophthora spp. in Cacao (Theobroma cacao L.). *G3 Genes/Genomes/Genetics* 2020;10:1713-25.

69. Karlström A, Gómez-Cortecero A, Nellist CF et al. Identification of novel genetic regions associated with resistance to European canker in apple. *BMC Plant Biol* 2022;22:452.

70. Tegtmeier R, Cobb-Smith D, Zhong G-Y et al. Identification and marker development of a moderate-effect fire blight resistance QTL in M. sieversii, the primary progenitor of domesticated apples. *Tree Genetics & Genomes* 2023;19:50.

71. Nutter FF. The role of plant disease epidemiology in developing successful integrated disease management programs. *General Concepts in Integrated Pest and Disease Management*. Springer, 2007, 45-79.

72. Djabali Y, Rincent R, Martin M-L et al. Plasticity QTLs specifically contribute to the genotype × water availability interaction in maize. *Theor Appl Genet* 2023;136:228.

73. Donoso JM, Picañol R, Serra O et al. Exploring almond genetic variability useful for peach improvement: mapping major genes and QTLs in two interspecific almond × peach populations. *Mol Breeding* 2016;36:16.

74. Marimon N, Luque J, Arús P et al. Fine mapping and identification of candidate genes for the peach powdery mildew resistance gene Vr3. *Horticulture Research* 2020;7:175.

75. Pascal T, Aberlenc R, Confolent C et al. Mapping of new resistance (Vr2, Rm1) and ornamental (Di2, pl) Mendelian trait loci in peach. *Euphytica* 2017;213:132.

76. Viruel MA, Madur D, Dirlewanger E et al. Mapping Quantitative Trait Loci Controlling Peach Leaf Curl Resistance. *IV International Peach Symposium* 1997;465:79-88.

77. Kobe B, Kajava AV. The leucine-rich repeat as a protein recognition motif. *Current Opinion in Structural Biology* 2001;11:725-32.

78. Dangl JL, Jones JDG. Plant pathogens and integrated defence responses to infection. *Nature* 2001;411:826-33.

79. Han G-Z. Origin and evolution of the plant immune system. *New Phytologist* 2019;222:70-83.

80. Gottin C, Dievart A, Summo M et al. A new comprehensive annotation of leucine-rich repeat-containing receptors in rice. *The Plant Journal* 2021;108:492-508.
81. Dodds PN, Rathjen JP. Plant immunity: towards an integrated view of plant-pathogen interactions. *Nat Rev Genet* 2010;**11**:539–48.

82. Sett S, Prasad A, Prasad M. Resistance genes on the verge of plant-virus interaction. *Trends in Plant Science* 2022;**27**:1242–52.

83. Rehder A. *Manual of Cultivated Trees and Shrubs*. The Macmillan Company, 1949.

84. Miller AJ, Gross BL. From forest to field: Perennial fruit crop domestication. *American Journal of Botany* 2011;**98**:26.

85. Doré C, Varoquaux F. *Histoire et amélioration de cinquante plantes cultivées*. Editions Quae, 2006.

86. Arús P, Gardiner S. Genomics For Improvement Of Rosaceae Temperate Tree Fruit. In: Varshney RK, Tuberosa R (eds.). *Genomics-Assisted Crop Improvement: Vol 2: Genomics Applications in Crops*. Dordrecht: Springer Netherlands, 2007, 357–97.

87. Aranzana MJ, Abbassi E-K, Howad W et al. Genetic variation, population structure and linkage disequilibrium in peach commercial varieties. *BMC Genet* 2010;**11**:1–11.

88. Brachi B, Morris GP, Borevitz JO. Genome-wide association studies in plants: the missing heritability is in the field. *Genome Biology* 2011;**12**:232.

89. Fetch T, McCallum B, Menzies J et al. Rust diseases in canada. *Prairie Soils and Crops* 2011;**4**:87–96.

90. Adaskaveg J, Soto-Estrada A, Förster H et al. *Peach Rust Caused by Tranzschelia Discolor in California*. UCANR Publications, 2000.

91. Parveaud CE, Broquaire J, Mercier V et al. Apricot susceptibility to blossom brown rot (Monilia spp.) and leaf rust (Tranzschelia spp.) under low-input production system. *15 international conference on organic fruit-Growing* 2020:427.

92. Giordani E, Nin S, Nencetti V et al. Inheritance for resistance to leaf curl [Taphrina deformans (Berk.) Tul.] in peach cultivars estimated by exploring mixed models. *Euphytica* 2022;**219**:1.

93. Iwata H, Minamikawa MF, Kajiya-Kanegae H et al. Genomics-assisted breeding in fruit trees. *Breeding Science* 2016;**66**:100–15.

94. Myles S, Peiffer J, Brown PJ et al. Association Mapping: Critical Considerations Shift from Genotyping to Experimental Design. *The Plant Cell* 2009;**21**:2194–202.

95. Bazakos C, Hanemian M, Trontin C et al. New Strategies and Tools in Quantitative Genetics: How to Go from the Phenotype to the Genotype. *Annual Review of Plant Biology* 2017;**68**:435–55.

96. Serrie M, Blanc A, Brun L et al. Adapting stone fruit trees to pesticide reduction: screening apricot and peach core collections for multi-pest resistance via integrative phenotyping under low pesticide management. *Acta Hortic* 2023:329–38.

97. Rife TW, Poland JA. Field Book: An Open-Source Application for Field Data Collection on Android. *Crop Science* 2014;**54**:1624–7.
98. Danecek P, Auton A, Abecasis G et al. The variant call format and VCFtools. *Bioinformatics* 2011; **27**:2156-8.

99. Hill WG, Weir BS. Variances and covariances of squared linkage disequilibria in finite populations. *Theoretical Population Biology* 1988; **33**:54-78.

100. Yu J, Pressoir G, Briggs WH et al. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. *Nature Genetics* 2006; **38**:203-8.

101. Laporte F, Charcosset A, Mary-Huard T. Efficient ReML inference in variance component mixed models using a Min-Max algorithm. *PLOS Computational Biology* 2022; **18**:e1009659.

102. Frichot E, Mathieu F, Trouillon T et al. Fast and Efficient Estimation of Individual Ancestry Coefficients. *Genetics* 2014; **196**:973-83.

103. Rincent R, Moreau L, Monod H et al. Recovering Power in Association Mapping Panels with Variable Levels of Linkage Disequilibrium. *Genetics* 2014; **197**:375-87.

104. Vitezica ZG, Legarra A, Toro MA et al. Orthogonal Estimates of Variances for Additive, Dominance, and Epistatic Effects in Populations. *Genetics* 2017; **206**:1297-307.

105. Roth M, Beugnot A, Mary-Huard T et al. Improving genomic predictions with inbreeding and nonadditive effects in two admixed maize hybrid populations in single and multienvironment contexts. *Yu J (ed.). Genetics* 2022; **220**:iyac018.

106. Li J, Ji L. Adjusting multiple testing in multilocus analyses using the eigenvalues of a correlation matrix. *Heredity* 2005; **95**:221-7.

107. Cinar O, Viechtbauer W. The poolr Package for Combining Independent and Dependent p Values. *Journal of Statistical Software* 2022; **101**:1-42.

108. Tsepilov YA, Shin S-Y, Soranzo N et al. Nonadditive Effects of Genes in Human Metabolomics. *Genetics* 2015; **200**:707-18.

109. van Eeuwijk FA, Bink MC, Chenu K et al. Detection and use of QTL for complex traits in multiple environments. *Current Opinion in Plant Biology* 2010; **13**:193-205.

110. Mangin B, Siberchicot A, Nicolas S et al. Novel measures of linkage disequilibrium that correct the bias due to population structure and relatedness. *Heredity* 2012; **108**:285-91.
Figures

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Phenotypic distribution and correlation between Max and AUDPC for A) peach rust and B) apricot rust. Distributions of Max are shown at the top of each graph, while the distributions of AUDPC are shown on the right.</figcaption></figure>

AV: Avignon, ET: Étoile-sur-Rhône, TO: Torreilles, SM: Saint-Marcel-lès-Valence and BE: Bellegarde.
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** A) Stacked bar plots with the variance of the random effects of genotype, genotype by location, genotype by year and genotype by year by location interactions, and residuals; calculated from the model presenting in Eq. (2) of Supplementary Method 3 and B) broad sense heritability calculated from Eq. (3) Supplementary Method 3.</figcaption></figure>

*Biotic stress observed only one year

**Biotic stress observed only in one location
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

**Figure 3.** GWAS results on main additive and non-additive genetic effect with A) QTL summary showing the physical position of GWAS significant SNPs. Each horizontal line contains QTLs of one biotic stress, organized by additive or non-additive analysis and either for the maximum of damage score (called Max) or the AUDPC. The dot colours are proportional to the -log₁₀(*P*-value) of the QTL top SNP. The x-axis indicates the physical positions on the peach or apricot genome. Boxplots illustrating the effect of a (B) dominant, (C) recessive and (D) overdominant SNP with a significant association for respectively shot hole (Max), blossom blight (Max) and apricot rust (AUDPC).
<figure><img src="image_8.png" /><figcaption>**Figure 4.** QTL summary showing the physical positions of GWAS significant SNPs for environment-specific models. Each horizontal line contains QTLs of one biotic stress, organized by environment and either for the maximum of damage score (called Max) or the AUDPC (except for blossom blight, where AUDPC was not available). The coloured dots are proportional to the -log₁₀(P-value) of the QTL top SNP. The x-axis indicates the physical positions on the genome (peach or apricot depending on trait).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 5.** QTL summary showing the physical position of GWAS significant SNPs for the two models dissecting the GxE interactions (MTMM and Meta-GWAS). Each horizontal line contains QTLs of one biotic stress, organized by GWAS model (Meta-GWAS and the 3 sub-model of MTMM: full, common and gei). The dot colours are proportional to the -log₁₀(*P*-value) of the QTL top SNP and the dot shape represent the considered variable (Max or AUDPC). The x-axis indicates the physical positions on the peach or apricot genome.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Fine characterization of two high confidence QTLs. A) Representation of the local LD between the top SNP of QTL 'Leaf_curl_4' and all other SNPs of the QTL. The X-axis represents the physical distance from the top SNP, while the Y-axis represents the r², corrected by the relatedness and structure, between SNPs and the top SNP. Each dot represents a SNP which is coloured according the p-value of Meta-GWAS. The vertical orange dot lines represent the QTL interval readjusted using the local LD. B) genes with relevant functions located in the QTL, with candidate genes in green and genes located in the QTL but not in local LD with the top SNP in red. C) Heatmaps of pairwise LD estimates within the genomic window around candidate gene 'PruarM.6G017900'. Top SNPs of the different high confidence QTLs are shown with red dot line and D) G-BLUPs distributions of apricot rust Max according to the haplotype constituted by the SNPs 'S6_974955','S6_975427','S6_975502' and 'S6_976066'. Statistical differences between groups were computed by pairwise Wilcoxon tests.</figcaption></figure>