# Insertion of a mMoshan transposable element in PpLMI1, is associated with the absence or globose phenotype of extrafloral nectaries in peach [Prunus persica (L.) Batsch

Patrick Lambert, Carole Confolent, Laure Heurtevin, Naïma Dlalah, Véronique Signoret, Bénédicte Quilot-Turion, Thierry Pascal

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

## To cite this version:

Patrick Lambert, Carole Confolent, Laure Heurtevin, Naïma Dlalah, Véronique Signoret, et al.. Insertion of a mMoshan transposable element in PpLMI1, is associated with the absence or globose phenotype of extrafloral nectaries in peach [Prunus persica (L.) Batsch. Horticulture research, 2022, 9, pp.1-13. 10.1093/hortre/uhab044. hal-03539501

HAL Id: hal-03539501

https://hal.inrae.fr/hal-03539501

Submitted on 23 Jun 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
**Title:** Insertion of a *mMoshan* transposable element in *PpLMI1*, is associated with the absence and globose phenotype of extrafloral nectaries in peach [*Prunus persica* (L.) Batsch]

**Running title:** a *mMoshan* is associated with the absence of EFNs in peach

**Authors :**

*Patrick Lambert¹*

Carole Confolent¹, ²

Laure Heurtevin¹

Naïma Dlalah¹

Véronique Signoret¹

Bénédicte Quilot-Turion¹

Thierry Pascal¹

**Affiliations**

¹INRAE, GAFL, Montfavet, F-84143, FRANCE

²INRAE, UMR GDEC, Clermont-Ferrand, F-63100, France

**Authors email adresses**

patrick.lambert.ugafl@inrae.fr

carole.confolent@inrae.fr

laure.heurtevin@inrae.fr

naima.dlalah@inrae.fr

veronique.signoret@inrae.fr

benedicte.quilot-turion@inrae.fr

thierry.pascal84@orange.fr (previous address: thierry.pascal@inrae.fr)

**Author for correspondence**

Patrick Lambert (Telephone +33 432 722 815; Fax +33 432722 702)
# 1 Abstract

Most commercial peach [Prunus persica (L.) Batsch] cultivars have leaves with extrafloral nectaries (EFNs). Breeders have selected this character over time, as they observed that the eglandular phenotype resulted in high susceptibility to peach powdery mildew, a major disease of peach trees. EFNs are controlled by a Mendelian locus (*E*), mapped on chromosome 7. However, the genetic factor underlying *E* was unknown. In order to address this point, we developed a mapping population of 833 individuals derived from the selfing of 'Malo Konare', a Bulgarian peach cultivar, heterozygous for the trait. This progeny was used to investigate the *E*-locus region, along with additional resources including peach genomic resequencing data, and 271 individuals from various origins used for validation. High-resolution mapping delimited a 40.6 kbp interval including the *E*-locus and four genes. Moreover, three double-recombinants allowed identifying *Prupe.7G121100*, a *LMI1-like* homeodomain leucine zipper (HD-Zip) transcription factor, as a likely candidate for the trait. By comparing peach genomic resequencing data from individuals with contrasted phenotypes, a MITE-like transposable element of the hAT superfamily (*mMoshan*) was identified in the third exon of *Prupe.7G121100*. It was associated with the absence and globose phenotype of EFNs. The insertion of the transposon was positively correlated with enhanced expression of *Prupe.7G121100*. Furthermore, a PCR marker designed from the sequence-variants, allowed to properly assign the phenotypes of all the individuals studied. These findings provide valuable information on the genetic control of a trait poorly known so far although selected for a long time in peach.

# 20 Introduction

Most peach [Prunus persica (L.) Batsch] cultivars have extrafloral nectaries (EFNs), or leaf-glands, on the leaf petioles, stipules, or margins¹,². EFNs are nectar-secreting glands, physically apart from the flowers. They have been observed on a vast diversity of species spanning over 93 families and 332 genera³,⁴. EFNs are mainly known for providing plants with indirect defense against herbivores and fungi, by attracting beneficial predatory arthropods, predominantly ants, and fungivorous mites with their sugary secretions³,⁵,⁶. EFNs can enhance plant-mite mutualisms by increasing mite abundance in
domatia, indirectly decreasing pathogen load. Therefore, predatory mites keep leaves free of microscopic herbivores while fungivorous mites clean the leaves of detrimental fungi⁷. EFNs may thus increase the potential of EFN-bearing peach cultivars to be protected from damaging organisms by naturally-occurring biological control agents⁵,⁶. From an extensive study of the main varieties of the peach, Gregory¹ observed that, for the great majority, gland shapes were well defined and that, in many cases, their shape could serve to separate groups of varieties. Indeed, gland shape was generally homogenous on typical shoots, although some cultivars could exhibit mixed glands. This author identified four main types of leaves, those with reniform (kidney-shape) glands, those with globose glands, glandless leaves and leaves having indistinctive glands. He reported that glands varied in number over the leaves of a same tree and were smaller on the leaf-margin than on the petiole. Gregory¹ also observed that reniform glands were associated with single crenate leaf-margins, whereas leaf-margins were doubly and deeply serrated in eglandular individuals. In the past, fruit breeding programs had inadvertently produced peach cultivars with glandless leaves, yet without determining the effects on either natural enemies or herbivorous pests²,⁸. Mathews et al.⁵ however, comparing glandular and eglandular peach trees derived from the selfing of the cultivar 'Lovell', observed that those trees with EFNs harbored significantly fewer herbivores than trees without EFNs. The latter also experienced lower growth and fruit production. Earlier, empirical observations showed that the absence of EFNs in peach cultivars resulted in high susceptibility to peach powdery mildew (PPM), one of the major diseases of the peach⁹,¹⁰. Additionally, in wild grape, Weber et al.⁷ demonstrated that adding foliar sugar to plant leaves increased the number of mutualistic mites inhabiting leaf domatia, and this was negatively correlated with the extent of the establishment of grape powdery mildew, a fungal disease similar to PPM. PPM is caused by *Podosphaera pannosa* var. *persicae*¹¹, a member of the Ascomycete fungi, which can be responsible for serious damages in peach orchards. Indeed, the disease may induce necrosis and malformation resulting in unmarketable fruits, premature drop and shoot stunting¹². For this reason, eglandular peach seedlings were systematically discarded during the selection process of most of the breeding programs. The Mendelian inheritance of the leaf-gland phenotype was first described by Connors¹³. The trait has an incomplete dominance, the absence being recessive and the globose shape of the nectaries representing the heterozygous phenotype. Further studies allowed to map the trait on a
single locus (*E*) on chromosome 7 of the peach<sup>14,15</sup>, but without identifying any factor responsible for the trait and its variations. Furthermore, this same region was found associated with a minor Quantitative Trait Loci (QTL) for resistance to PPM in *P. ferganensis*<sup>14,16</sup>. These various studies contributed to provide evidence that EFNs might play a role in lowering PPM incidence and could be of most interest for limiting the populations of some classes of detrimental herbivores and fungi in the trees. Therefore, further investigations deserved to be conducted to identify the factor underlying the *E* locus. For this reason, in-depth study of the *E* locus was carried out, in the frame of our breeding program for resistance to pests and diseases in peach. The main objectives of the current work were to develop a high-resolution map of the *E* locus, then investigate the underlying genomic region in order to identify the factor involved in the variation of the leaf-gland phenotype as well as its possible link to susceptibility to PPM, apart from the indirect defense to fungi provided by EFNs. Then, accessorily, develop PCR marker(s) to facilitate early selection of glandular seedlings. With this aim, a large mapping population of 833 individuals, referred to as 5392<sup>2</sup>, was developed from the selfing of 'Malo Konare' (clone S5392), a canning peach cultivar with globose leaf-glands, from Bulgarian origin<sup>17</sup>. This cultivar was selected as it was heterozygous for the trait and part of our breeding program for resistance. Peach genomic resequencing data from contrasted cultivars were used for in-depth investigation of the *E*-locus region.

Additional resources including offspring derived from another cross, as well as a collection of contrasted cultivars from various origins were used to support our findings. The outcomes of this study will provide valuable information on a trait little studied in peach so far and more widely in *Prunus* species.

Furthermore, they would benefit our breeding program aimed at developing multi-resistant elite peach cultivars.

## 76 Results

### 77 Phenotypic evaluation

Seven hundred and seventy-nine progenies out of the initial 833 of the 5392² were observed over two years, among which 197 from the initial population. Two hundred and two (26%) were eglandular, 382 (49%) globose and 195 (25%) reniform. This distribution was in agreement with the (1:2:1) segregation ratio expected for a Mendelian trait in this type of population (χ² =0.41). Regarding the BC2, 62
individuals had globose leaf-glands and 60 had reniform leaf-glands. As regards the collection, 112 cultivars and the two wild peach relatives were scored reniform, 30 globose, four eglandular and one indeterminate (Table S1). For *Prunus kansuensis* S1429, a few phenotypic differences with the other accessions were observed: leaves were homogenous but EFNs were included in the margin of the lower part of the leaf-blade instead of the upper ridge of the petiole. Regarding PER2.3N#1 (S7314), a possible triploid scored indeterminate, no regular leaf-gland was noticeable but a number of small picks on the petiole, close to the leaf-blade. Finally, with respect to leaf-margins, a close association was observed between deeply serrated leaves and the eglandular phenotype in the population 5392². Eglandular individuals had sharp doubly well-defined leaf serrations contrary to those with globose or reniform glands, which had leaves with rounded, shorter crenellations. Crenellations were generally slightly more pronounced in globose individuals (Fig. 1). Regarding the collection, the same association was observed for the four eglandular accessions as compared to the others.

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Photographs of the three types of leaves observed in the population 5392². Top photos show the three different phenotypes observed for the EFNs. Bottom photos show the leaf-margins associated with each of the above phenotypes. (A) Eglandular 'S10215', (B) Globose 'Malo Konare' (C), Reniform 'S10216'.</figcaption></figure>
99

100 Genetic map of linkage group 7 and high-resolution mapping of the *E* locus

101 The map of G7 derived from the 212 initial individuals covered a total genetic distance of 80.9 cM (Fig.

102 2) spanning a physical distance of 19,892,186 bp (88.85% of chromosome 7).

103

G7_Malo_Konare

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

104

105 Fig. 2 Genetic map of linkage group 7 of ‘Malo Konare’ developed from the initial mapping
106 population of 212 individuals. The EFN locus (E) is in bold and in italics. Genetic distances are in
107 centiMorgan (cM).

108

109 The map was composed of 18 SNPs among which six, including ASPP900, collocated with the *E* locus,
110 at 49.6 cM, spanning a physical distance of 107.8 kbp. The physical distance between the SNPs on either
111 side of the *E*-locus region (SNP_IGA_776067 and SNP_IGA_777469) was 665.7 kbp for a genetic
112 distance of 1.6 cM. No significant deviation of marker segregation was observed (*P* < 0.05). In addition
131 Table 2 Predicted genes observed in the 70.5-kbp genomic region comprised between

132 SNP_IGA_776067 and SNP_IGA_776214

<table><thead><tr><td>Annotated gene</td><td>Position on Peach v2.0</td><td>Swissprot description/match</td><td>TAIR description/match</td></tr></thead><tbody><tr><td>Prupe.7G120700</td><td>Pp07:14410619..14416964</td><td>Pyruvate kinase, cytosolic isozyme (Glycine max) /Q42806</td><td>Pyruvate kinase family protein/ AT3G52990.1</td></tr><tr><td>Prupe.7G120800</td><td>Pp07:14417865..14420155</td><td>Uncharacterized</td><td>Sequence-specific DNA binding transcription factors/ AT3G10040.1</td></tr><tr><td>Prupe.7G120900</td><td>Pp07:14425004..144427458</td><td>Uncharacterized</td><td>Glycoside hydrolase family 28 protein/ AT2G33160.1</td></tr><tr><td>Prupe.7G121000</td><td>Pp07:14428432..14431463</td><td>F-box protein PP2-A15 (Arabidopsis thaliana) /Q9LF92</td><td>Phloem protein 2-A15/ AT3G53000.1</td></tr><tr><td>Prupe.7G121100</td><td>Pp07:14436305..14437630</td><td>Putative homeobox-leucine zipper protein ATHB-51 (Arabidopsis thaliana) /Q9LZR0</td><td>Homeobox 51/ AT5G03790.1</td></tr><tr><td>Prupe.7G121200</td><td>Pp07:14438623..14440853</td><td>60S ribosomal protein L24 (Prunus avium) /Q9FUL4</td><td>Ribosomal protein L24e family protein/ AT3G53020.1</td></tr><tr><td>Prupe.7G121300</td><td>Pp07:14441874..14445148</td><td>Protein kinase dsk1 (Schizosaccharomyces pombe)/P36616</td><td>Ser/arg-rich protein kinase 4/ AT3G53030.1</td></tr><tr><td>Prupe.7G121400</td><td>Pp07:144459229..14449576</td><td>Uncharacterized</td><td>Sulfite exporter TauE/SafE family protein 4/ AT2G36630.1</td></tr><tr><td>Prupe.7G121500</td><td>Pp07:14467931..14469921</td><td>Embryonic protein DC-8 (Daucus carota)/P20075</td><td>Embryonic cell protein 63/ AT2G36640.1</td></tr><tr><td>Prupe.7G121600</td><td>Pp07:14471118..14474276</td><td>Probable glycosyltransferase At5g03795 (Arabidopsis thaliana) /Q9FFN2</td><td>Exostosin family protein/ AT5G03795.1</td></tr><tr><td>Prupe.7G121700</td><td>Pp07:14480000..14482717</td><td>Pentatricopeptide repeat-containing protein At5g03800 (Arabidopsis thaliana) /Q9FFN1</td><td>Pentatricopeptide repeat (PPR) superfamily protein/ AT5G03800.1</td></tr><tr><td>Prupe.7G121800</td><td>Pp07:14483790..14494599</td><td>E3 ubiquitin-protein ligase UPL7 (Arabidopsis thaliana) /Q9SCQ2</td><td>Ubiquitin-protein ligase 7/ AT3G53090.2</td></tr></tbody></table>
Reads from the eglandular 'S10215', the reniform 'S10216', 'Summergrand' and 'Pamirskij 'S10215', as well as the globose 'Zephyr' and 'Malo Konare' were aligned onto Peach v2.0.a1 derived from the reniform peach cv. Lovell (Plov2-2N) and compared. A total of two hundred and seventy-seven variants between 'S10215' and 'S10216'and heterozygous in 'Malo Konare', were identified among which six SNPs and an indel in the ASPP899-ASPP901 region (Table S2). However, no relationship was observed between any of the variants and the trait, except for the indel, which clearly differentiated eglandular, reniform and globose accessions. For the other 276 variants, the eglandular 'S10215' had the same haplotype as the reniform 'Summergrand' and 'Lovell' (Plov2-2N), as well as the globose 'Zephyr'. In contrast, the reniform 'S10216' was highly similar to 'Pamirskij 5' (reniform), except for an 11-kbp region upstream of the indel, for which 'Pamirskij 5' had the same haplotype as the above four other accessions (Table S2). The indel was located in *Prupe.7G121100*, a gene annotated as putative homeobox-leucine zipper protein ATHB-51 (Table 2). According to Gene Ontology, *Prupe.7G121100*, has a DNA-binding transcription factor activity and is involved in bract formation and leaf morphogenesis. Based on these findings, 25 primers (Table S3) were developed from consensus regions between 'S10215' and 'S10216' in order to sequence the interval encompassing *Prupe.7G121100*, as well as the 100-bp gap which remained in the peach genome sequence reference⁴⁴ (Peach v2.0.a1) immediately upstream of the CG (position Pp07:14436205..14436304). The sequencing of the gap region resulted in sequences 9-fold longer than expected (905 bp and 903 bp for S10215 and S10216 respectively), therefore impacting coordinates downstream (Fig. S1). Sequence comparison allowed identifying a 590-bp insertion in the last coding DNA sequence (CDS) of *Prupe.7G121100*, in the eglandular 'S10215', as well as two additional polymorphisms due to differences in the number of CT repeats in two SSRs present in the gap region (Fig. S1). The 590-bp insertion was located between positions Pp07:14437331 and Pp07:14437332, disrupting the initial reading frame (Fig. S1). BLASTN search against NCBI database allowed finding a high similar hit (98% of identity) with an insertion fragment of 588
160 bp, upstream of the start codon of a chalcone isomerase (CHI) gene of peach (Sequence ID: KF990613.1). This insertion was identified as a MITE-like *Moshan* (*mMoshan*) transposable element of the hAT superfamily. BLASTN search with the sequence inserted in *Prupe.7G121100*, against Peach v2.0.a1, returned 91 additional highly-similar hits (> 95% identity) spanning all the chromosomes, all starting from the third 5' nucleotide of the inserted element. The most similar (100% of identity from the third 5' nucleotide) was located on chromosome 5 (Pp05:16,628,569-16,629,156), 460 bp upstream of the start codon of *Prupe.5G208500*, a homolog of AGL8 (agamous-like 8) transcription factor. Nevertheless, a fine analysis of the insertion sequence highlighted some differences with the other transposable elements. The 92 above transposons had terminal inverted repeats (TIRs) composed of 13 or 14 complementary nucleotides and 8-nucleotide target site duplication (TSD). Regarding the insertion, the TIR in 3' was composed of 13 nucleotides identical to that of 90 of the 92 transposons. However, only 10 nucleotides of the 5' TIR were complementary with those of the 3' TIR (Fig. S1). In addition, no direct repeat sequence was observed at the target insertion site and therefore no TSD. A likely hypothesis is a deletion in the original sequence (GACGAGCCTAGGGGTGGGCAC) where "GACGAGCC" was the TSD, the deleted motif "CGAGCCTAGG" and the original 5' TIR started with the motif "TAGGG".

177

178 Analysis with FGENESH

179 The analysis with FGENESH was performed for both variants, using the genomic sequence of *Prupe.7G121100* supplemented by the sequence of the gap region, and various dicot plant species as models, among which *P. persica*, *A. thaliana*, *M. domestica* and *L. esculentum*. One single prediction was obtained with the reniform sequence (Fig. S2). The addition of the gap-region sequence leaded to a primary transcript composed of three CDS instead of two, the initial start codon being replaced by another one 431 bp upstream. Queries of GDR_*RefTransV1* and NCBI database using the
185

sequence of the resulted transcript, validated the prediction, three transcripts (P.persica_gdr_reftransV1_0044698, *P. dulcis* LOC117635596 and *P. avium* LOC 110754228) having sequences highly similar with the predicted transcript (100%, 99% and 98% identity respectively). For eglandular accessions, four different predictions were obtained, all of these including an additional CDS in the 3' region. Differences between predictions were linked to the proportion of the transposon included in the third CDS and the position of the fourth.

186

187

188

189

190

191

## Expression analysis of *Prupe.7G121100*

192

193

194

Relative expression levels of *Prupe.7G121100* in leaves were assessed in three eglandular, two globose and three reniform cultivars, as well as two wild species, *P. davidiana* P1908 and *Prunus kansuensis* S1429 (Fig. 4).

195

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

196

197

198

199

200

201

**Fig. 4 Relative expression of *Prupe.7G121100* in ten accessions contrasting for EFNs.** Expressions were normalized with the constitutive genes *PpTEF2* and *PpRPL13*. Eglandular, globose and reniform individuals are denoted by the letter E, G or R, before the accession number, respectively. The three eglandular individuals are framed red. The expression of *P. kansuensis* is represented by two samples: R-S1429M (leaf-margin) and R-S1429P (upper-petiole region).

202
Contrary to our initial expectations, *Prupe.7G121100* had significantly higher expression levels (p < 0.001) in eglandular accessions, than in both reniform and globose individuals, either before or after normalization with *PpTEF2* and *PpRPL13*. Normalized differential expressions were comprised between 0.049 ± 0.0052 (mean ± SE) and 0.1454 ± 0.03675 for reniform individuals, 0.1423 ± 0.01 and 0.1778 ± 0.0093 for the globose ones, and between 0.6846 ± 0.0394 and 1.1818 ± 0.0627 for eglandular accessions (Fig. 4), thus showing a negative correlation between the expression level of *Prupe 7G121100* (p < 0.001) and the presence of EFNs. No significant difference was observed between the two samples of *Prunus kansuensis* S1429 (p < 0.01). In comparison, differential expression values before normalization, were comprised between 1 ± 0.10 (mean ± SE) and 2.96 ± 0.75 for reniform individuals, 2.89 ± 0.19 and 3.62 ± 0.20 for the globose ones, and between 13.93 ± 0.77 and 23.79 ± 0.99 for eglandular accessions.

## 3' RACE PCR and comparison of the alleles of the transcript

Nested PCRs based on sense primers associated with the AUAP antisense primer gave single amplicons for the reniform accessions only, whereas those carried out on eglandular accessions produced mixtures of amplicons of different sizes (smears). In contrast, those carried out using the antisense primers developed from each of the four predictions, (Table S6) gave the expected results, with amplicons present or absent according to the prediction considered. Sequences derived from the amplicons confirm the insertion of the 179 first nucleotides of the transposon after position 134 of the third exon of the initial transcript, as well as the presence of a fourth exon in the eglandular accessions (Fig. S1 and Fig.3).
This confirms that prediction #3, which includes *P. persica* in the model, is the only valid (Fig. S2).
Regarding the reniform accessions, the transcript was as expected. The size of the eglandular transcript was 99 nucleotides longer than that of the reniform one (777 and 678 nucleotides respectively) resulting in a larger predicted protein (258 and 225 amino acids respectively). Moreover, major changes were observed in the eglandular transcript : 33 amino acids of the 3' end were replaced by 65 others in the eglandular transcript (Fig. S2).
<figure><img src="image_6.png" /><figcaption>**Fig. 3 Diagrams of PpLMI1.** Spliced transcripts are displayed below their respective primary transcripts (A) Truncated *PpLMI1* (Prupe 7.G121100) as annotated in Peach v2.0.a1. CDSs are shown as blue rectangles, 5' and 3'-UTR as white rectangles, introns as black lines. The start and stop codons are shown as yellow and green rectangles respectively. Upper coordinates represent current positions on Peach v2.0.a1 (Pp07) although they are no longer relevant as the gap upstream of Prupe 7.G121100 is longer than indicated. Lower coordinates represent the distance from the first nucleotide of the 5'-UTR. (B) *PpLMI1* as observed in reniform individuals. (C) *PpLMI1* as observed in eglandular individuals. Regions not shared with the reniform transcript are shown in red and correspond to transposon segments (coding sequences, intron and 5'-UTR).</figcaption></figure>

## Genotyping with the ASPP900 marker

One thousand and fifty individuals in total were genotyped with the ASPP900 marker, among which two hundred and seventy-one individuals used for validation, including 149 accessions (Table S1). For all of them except 'S7314', genotypes were consistent with phenotypes and globose individuals could also be clearly differentiated from reniform ones. 'S7314' was considered as a possible triploid derived from the eglandular 'Prosser 2.1N'; however, it was genotyped as globose (heterozygous) and phenotyped as undifferentiated. These discrepancies do not question the efficiency of ASPP900, but are
247

rather due to its peculiar genotype. In addition, this raises doubts on the single parental origin of this
accession.

248

249

Discussion

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

The aim of our study was to identify the genomic factor responsible for the presence/absence of EFNs in peach, a Mendelian trait previously mapped on chromosome 7 (ref. 14, 15), but little studied so far.

The fine-mapping approach allowed delimiting the trait to an interval of 10.7 kbp between positions Pp07:14,428,469 and Pp07:14,439,211. Micheletti et al.¹⁸, using the ISPC 9K SNP peach array¹⁹ and a collection comprising 750 reniform and 190 globose accessions, identified a single SNP, SNP_IGA_776161 (position Pp07:14,469,094) as associated to the leaf-gland type. This association was not fully congruent with our observations as ‘S10215’ (eglandular), ‘Zephyr’ (globose), ‘Summergrand’, ‘Rubira’ and the peach genome reference derived from ‘PLov2-2N’ (reniform) had the same allele combination (C/C) for this SNP, whereas ‘S10216’ and ‘Pamirskij 5’, both reniform were T/T. Nevertheless, taking into account the limited number of SNPs on the array corresponding to the trait region, as well as possible misclassification of some individuals of the collection, the results of the two studies were convergent. Coupling the results of the fine-mapping approach with the comparison of the genomic sequences of accessions contrasting for the trait, then allowed to clearly identifying a single candidate gene, *Prupe.7G121100*, among the three genes included in the above interval, and more broadly, among the 12 genes comprised in the longer 70.5-kbp genomic region encompassing the latter.

Indeed, *Prupe.7G121100* has two variants: the regular one, associated with the presence of EFNs and homozygous in reniform individuals, and a second one including a 590-bp insertion homozygous in eglandular individuals. This insertion was identified as a MITE-like transposable element of the hAT superfamily, termed as *mMoshan*²⁰. Miniature inverted-repeat transposable elements (MITEs) are non-autonomous class II transposable elements. They are considered a major driving force for generating allelic diversity in plant genomes²¹. MITES account for 3.89% of the peach genome²² with 0.16% for the 491 *Moshan* elements identified. *Moshan* elements are unique to Rosaceae and the *mMoshan* class is predominant with 432 elements²⁰. Interestingly, two of these *mMoshan* elements generated no obvious target site duplication, as the element inserted in *Prupe.7G121100*, suggesting that these three elements
were atypical. Wang et al.²⁰ identified 29 *mMoshan* which were inserted in genes, among which 14 in exons. The 29 genes were distributed over all the chromosomes but none in chromosome 7. The *mMoshan* in *Prupe.7G121100* was not detected probably because Peach genome v2.0.a1 was derived from the reniform double haploid 'Lovell' Plov2-2N and therefore does not include the inserted element.
This author observed that genes including *mMoshan* elements showed relatively lower expression levels compared with genes lacking these elements and this was consistent with previous studies on MITES²³.
However, this was not the case in our study since *Prupe.7G121100* demonstrated enhanced expression in eglandular individuals compared to that in globose and reniform ones, which were quite similar.
*mMoshan* elements contain several *cis*-regulatory elements such as MYB and WRKY binding sites in the first third of the sequence, which could be involved in upregulation of the transcription²⁰. However, when we take into account the minor differences in gene expression observed between reniform and globose individuals, and the similarity of the phenotypes of their leaf-margins, as well as the incomplete dominance of the trait, this seems not correlated. It would be therefore interesting to investigate possible functional differences between the two alleles This point needs a dedicated approach.
*Prupe.7G121100* was annotated as putative homeobox-leucine zipper protein ATHB-51, a member of the class I (HD-Zip I) superfamily of transcription factors. (HD-Zip) proteins are unique to plants. They include the peculiar combination of a DNA-binding homeodomain (HD) and an adjacent Leucine zipper (Zip) motif, which mediates protein-dimer formation²⁴. Saddic et al.²⁵ identified ATHB-51 as a meristem identity regulator and named it LATE MERISTEM IDENTITY (*LMII*) based on its regulation functions; accordingly, we will further refer to *Prupe.7G121100* as *PpLMII*. These authors showed that *LMII* was a direct target of *LEAFY* (*LFY*), a central meristem identity regulator in *Arabidopsis thaliana* as well as a direct upstream activator of a second meristem identity regulator, the MADS-box transcription factor *CAULIFLOWER* (*CAL*). *LMII* acts together with *LFY* to induce *CAL* expression, the interaction between these three genes corresponding to a feed-forward loop transcriptional network motif²⁶. *LMII* thus belongs to the complex of genes including others transcription factors, such as APETALA1 (*API*), involved in the meristem identity switch leading to flower formation²⁵,²⁷.
Interestingly, the *mMoshan* transposable element identified on chromosome 5, with the highest
percentage of identity with that inserted in *PpLMI1*, was in the promoter region of an homolog of AGL8 (agamous-like 8) transcription factor, a MADS-box negatively regulated by APETALA1, suggesting a possible involvement of APETALA1 in the regulation of *PpLMI1*. However, *LMI1* has also additional *LFY*-independent roles in leaf morphogenesis and bract formation²⁵. For instance, *LMI1* regulates leaf growth in *Arabidopsis thaliana*²⁸ as well as organ proportions such as stipules size, via an endoreduplication-dependent trade-off, that limits tissue size and cell proliferation, through the activation of the mitosis blocker *WEE1*²⁹. Moreover, modifications of *GhLMI1-D1b*, one of its homologues³⁰ were found to be responsible for the major leaf shapes in Upland cotton (*Gossypium hirsutum*). This designates *LMI1-LIKE* genes (along with the *KNOXI* genes) as evolutionary hotspots that have been recruited in angiosperms to modify leaf shape³¹. In this way, Chang et al. ³² demonstrated that *LMI1-like* and *KNOX1* genes coordinately control leaf development in dicotyledons and that different expression patterns of these two genes correspond to the formation of different leaf marginal structures. The same way, loss of function of *CrLMI1*, a likely ortholog of *LMI1* was reported to decrease leaf serration in *Capsella rubella*³³. This is in agreement with the results of our study, in which increased leaf-margin serration was found strictly associated with the absence of EFNs, concurrently with the higher expression of *PpLMI1*. This relationship between leaf serration and absence of EFNs was already reported in previous studies¹³. These findings thus contribute to confirm the possible involvement of *PpLMI1* in leaf-margin structures in peach and accordingly in the phenotype of EFNs. Likewise, in cucumber (*Cucumis sativus*), *mict*, a class I HD-Zip factor which sequence had 52% of identity with *LMI1*, regulates multicellular trichome development³⁴. Regarding the EFNs, however, molecular genetic understanding of their formation is still underdeveloped. No study to date is available in tree species and only a few ones have been published in annual plants. For instance, Hu et al.³⁵ identified *GaNEC1*, a gene encoding a PB1 domain-containing protein, as positive regulator of nectary formation in cotton, which silencing led to a smaller size of foliar nectary phenotype. However, EFNs were located in the leaf midribs and their conformation was different than peach EFNs. Phenotypic diversity usually results from diversity in the genetic organization, regulation and/or expression of underlying developmental programs⁴. In the case of EFNs, such underlying programs have been poorly investigated. The gene *CRABS CLAW (CRC)*, a YABBY transcription factor³⁶,³⁷, appears to be an early-
functioning regulator of the development of both floral and extrafloral nectaries in core eudicots³⁸,³⁹.

But while the location of floral nectaries may be determined by *CRC* along with several upstream MADS box floral homeotic genes and other unknown regulatory genes³⁸, the development of EFNs may involve the recruitment of different transcriptional control networks than those needed in floral nectaries³⁹. This means that the program involved in EFN development may be closely associated with that of the EFN-bearing organ⁴, the leaf, in the case of peach trees. As a result, the functional characteristics of *LMI1*, its involvement in leaf morphogenesis as well as in the meristem identity switch leading to flower formation, suggest that this transcription factor might have a pivotal role in the regulation of different characteristics of the leaf. This makes *PpLMI1* a most likely candidate for the presence/absence of EFNs in peach. Therefore, a plausible hypothesis is that functional modification of *PpLMI1* associated with the insertion of the HD-Zip I element might trigger endoreduplication. This would result in changes in the cell-wall composition of the lamina as well as in leaf margins, through a developmental program involving target-genes of *PpLMI1*, leading notably to serrated leaves and the absence of EFN. In addition, changes in cell-wall composition of the leaf-blade surface could thus make easier the development of fungi, such as *Podosphaera pannosa*, on the leaves. Modification of the cell walls at regions targeted by pathogen attack is a common response to infection and the inability to do so, or the presence of weakened cell walls, might explain, at least in part, the susceptibility to pathogens⁴⁰. As a result, these changes, along with the absence of the positive effects associated with the presence of domitia-inhabiting mutualist mites, and fungivore mites attracted by EFN nectar, might be responsible for enhanced susceptibility to PPM in eglandular individuals, as compared to those with EFNs. Further studies need however to be undertaken in order to assess our hypothesis.

350

## Conclusion

In this study, we were interested in identifying the genetic factor responsible for the presence/absence of EFNs in peach. In our knowledge, this is the first time that a molecular genetic approach has been undertaken to clarify the genetic basis of this Mendelian trait, in peach and, more broadly in *Rosaceae* perennial crops. Based on our results, *PpMLI1* appears the most likely candidate gene for this character.

A comprehensive study of the genomic region including *PpMLI1* study did not bring to light another
356

alternative candidate. In addition, its characteristics, regulation functions as a meristem identity
regulator as well as its role in leaf morphogenesis, make it highly plausible its involvement in the control
of the presence/absence of EFNs as well as its association, at some extent, with the variation of the
susceptibility to PPM, in link with cell-wall changes. However, this has to be further validated
functionally. Virus induced gene silencing (VIGS) method could be considered as a relevant approach,
as genetic transformation in peach is currently an obstacle. In addition, a broader study including the
expression of *PpMLI1* in the meristem as well as that of genes interacting with *PpMLI1* or target genes
such as *WEE1*, may be undertaken to elucidate the molecular interactions underlying this interesting
trait. This study is thus a first step. Nevertheless, in the short term and from a breeder perspective,
ASPP900 marker already allows differentiating the different phenotypes at the seedling level, and could
then be used in peach breeding programs.

357

358

359

360

361

362

363

364

365

366

367

## Material and methods

368

### Plant material

369

The initial mapping population included 212 individuals derived from the self-pollination of 'Malo
Konare' (clone S5392). 'Malo Konare' is a canning peach cultivar developed in 1984 at the Fruit-
growing Institute in Plovdiv (Bulgaria). It originated from the cross 'Stoika' × 'New Jersey Cling 97',
has globose leaf-glands and shows strong resistance to powdery mildew. 'Stoika' was for its part derived
from 'House Kling' and 'Ferganskyi Zheltyi' (1973), a clone of *Prunus ferganensis*. The population
was further extended to 833 individuals for the fine mapping of the leaf-gland region and identifying
recombinants. This population will be referred to as 5392². In addition, 271 individuals were used to
validate phenotype/genotype association in different genetic backgrounds: at first, 149 accessions with
contrasting leaf-gland phenotypes (Table S1), including 143 peach cultivars from various origins, two
accessions of wild species close to peach, *Prunus davidiana* (Carr.) and *Prunus kansuensis* (Koehne),
one accession of *Prunus ferganensis* (Kost. & Rjab.), two double haploids and a possible triploid; then
a sample of 122 individuals from a complex breeding population, referred to as BC2⁴¹. The latter was
derived from two successive crosses (F₁ and back-cross) including *Prunus davidiana* clone P1908 and
peach cv. 'Summergrand', followed by a final cross derived from a mixture of pollen of the back-cross

370

371

372

373

374

375

376

377

378

379

380

381

382
population and 'Zephyr' as maternal parent. These 271 individuals were planted in triplicate and grown in three different places: greenhouse and tunnels for the cultivars, orchards and tunnels for the BC2. All the individuals were conserved at the *Prunus* Biological Resource Center of INRAE in Montfavet, except the two double haploids and the possible triploid that were conserved at the *Prunus-Juglans* Biological Resource Center, Domaine des Jarres, 33210 Toulenne.

## 388 DNA isolation

Samples of young leaves from each of the individuals were collected in the spring. Genomic DNA was subsequently isolated using the Qiagen DNeasy 96 Kit (https://www.qiagen.com) according to the manufacturer's instructions. DNA of each sample was at first assessed for quality using a NanoDrop™ ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA USA) and then quantified using Quant-iT™ Picogreen® reagent (Invitrogen Ltd.2, Paisley UK). Stock solutions of genomic DNAs were then diluted to a final concentration of 40 ng/µl.

## 395 Leaf-gland phenotyping

Leaf-glands were observed over two years, on five to ten leaves from different parts of each of the trees (progenies and cultivars). Individuals were classified under the three phenotypes encountered: reniform, globose and eglandular (no leaf-gland observed). Those trees that were planted in triplicate were scored individually. Leaf-margins were examined concurrently to EFNs, as an association between the eglandular phenotype and deep leaf-serration was previously reported.

## 401 Next Generation Sequencing of accessions

Additionally to the reniform double-haploid peach reference 'Lovell' (PLov2-2N), which sequence is available at the GDR (https://www.rosaceae.org/species/prunus_persica/genome_v2.0.a1), seven peach accessions were used for genome comparison of the leaf-gland region: 'Summergrand', 'Pamirskij 5' and 'Rubira' (reniform), 'Zephyr' and 'Malo Konare' (globose) and two individuals derived from the self-pollination of 'Malo Konare', 5392<sup>2_60</sup> (eglandular) and 5392<sup>2_76</sup> (reniform) renamed 'S10215' and 'S10216' respectively. These seven accessions were sequenced by MGX GenomiX (Montpellier, France, http://www.mgx.cnrs.fr). In brief DNA libraries were prepared using the Nextera DNA Flex
409

Library preparation kit from Illumina (Illumina Inc. San Diego CA, USA) following recommendations
provided by the supplier. 125-bp paired-end sequencing was performed utilizing the Illumina HiSeq
2500 sequencing platform and the sequence by synthesis (SBS) technique. Base calling was performed
by the Real Time Analysis (RTA) software. Raw Illumina paired-end reads were subsequently trimmed
using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Potential contaminants
were investigated using FastQ Screen software (Babraham Institute) and Bowtie2 aligner (http://bowtie-
bio.sourceforge.net/bowtie2/index.shtml). Resulting reads were aligned onto Peach v2.0.a1 using
BWA-MEM (v0.7.12-r1039) and the Ppersica_298_v2.0.fa version. BAM files (*.sorted.bam and
*.sorted.bam.bai) were generated in order to visualize sequences under the Integrative Genomics Viewer
(IGV) tool⁴².

410

411

412

413

414

415

416

417

418

419

## Marker development and genotyping

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

'Malo Konare' has been genotyped earlier in the frame of the European project Fruitbreedomics¹⁸, using
the IPSC peach 9K SNP array v1¹⁹. Based on the available SNP dataset, a first set of heterozygous SNPs
was selected to develop the genetic map of linkage group 7 of 'Malo Konare' and insure a sufficient
coverage of the group. Genotyping was done using the PCR-based KASP™ (Kompetitive Allele
Specific PCR) method from LGC Biosearch technologies (https://www.biosearchtech.com/). Primer-
triplets (two competitive allele-specific forward primers and one common reverse primer for each
marker) were developed from the 60-bp genomic sequence available on either side of the SNPs
(https://www.rosaceae.org/species/rosaceae_family_genera/IRSC_SNP_array), using Primer3⁴³ under
the following primer-picking conditions: optimal size of the amplicons 75 bp (min 62 bp, max 85 bp),
Tm 65°C (min 55°C, max 72°C), primer size 25 bp (min 20 bp, max 32 bp), max self-complementarity
7, max 3' self-complementarity 3, left primer end 61 bp. Primers triplets were compared with Peach
v2.0.a1⁴⁴, using the Basic Local Alignment Search Tool (BLAST®) at the Genome Database for
Rosaceae (GDR: https://www.rosaceae.org/blast/). Those aligning to single positions were selected for
genotyping the starting mapping population (Table S4). In a second step, additional SNP markers
focused on the interval encompassing the *E* locus were developed in order to identify recombinant
individuals. This was done using Next-Generation Sequencing (NGS) data derived from 'Malo Konare'.
436

BAM files were aligned onto Peach v2.0.a1 and visualized with IGV⁴². The region containing *E* was
examined for SNP/indel discovery and reads including heterozygous SNP/indels compatible with the
KASP™ method were retrieved. Primer-triplets were then developed as above (Table S4). Mix
preparation and PCR reactions were performed using the KASP™ genotyping chemistry and conditions.

437

438

439

440

## Genetic map of linkage group 7

441

442

443

444

445

446

447

448

449

450

In a first stage, linkage group 7 (G7) of 'Malo Konare' was constructed using the mapping dataset derived from the SNP-set selected from Micheletti et al. ¹⁸. Genotypic data were coded as F₂-progeny type according to the JoinMap coding system. The leaf-gland trait was similarly coded as a co-dominant Mendelian trait. Linkage analyses were performed using JoinMap 4.1⁴⁵. The recombination fraction value was set at 0.4 and grouping was performed using the independence logarithm of odds (LOD) calculation function and a minimum LOD score threshold of 3. The Kosambi mapping function⁴⁶ was used to translate recombination frequencies into genetic distances. Linkage group 7 was established using regression mapping procedure with three rounds per sample. In a second stage, SNP markers developed for the high-resolution mapping in the interval including the *E*-locus region were added to the genotypic data file and mapped similarly.

451

## High-resolution mapping of the *E* locus

452

453

454

455

The extended population was genotyped using the SNP markers flanking the *E* locus in the genetic map.
Recombinant individuals in the interval were identified and genotyped with newly developed markers.
Individuals identified as recombinants in the new interval were genotyped again with a new marker-set.
This process was repeated iteratively until no further recombinant was observed.

456

## *In silico* analysis of the region encompassing the *E* locus

457

458

459

460

461

The genomic region delimited by the SNP-pairs, which allowed identifying the most informative recombinants, was analyzed for variants. This was done by aligning the sorted.bam files of 'Malo Konare', 'S10215' and 'S10216' onto Peach v2.0.a1 (Ppersica_298_v2.0.fa version), under IGV⁴², and by comparing them. Differences observed were then compared with sorted.bam files of 'Zephyr', 'Summergrand', 'Pamirskij 5' and 'Rubira' in order to check consistency of differences regarding leaf-
462 gland phenotype, in different genetic backgrounds. The genomic region defined above was examined
463 for the presence of predicted genes, using JBrowse on the Genome Database for Rosaceae
464 (https://www.rosaceae.org/jbrowse/). Positions of the observed differences were compared with those
465 of the genes and their sub-features, then, genomic sequences of the candidate genes (CGs), associated
466 transcripts and predicted protein sequences, homologies and gene functions were downloaded
467 (https://www.rosaceae.org/node/4017147). NGS reads corresponding to the position of the selected CG-
468 variants were retrieved for 'S10215' and 'S10216' using IGV⁴², imported into CLC Main Workbench
469 version 12 (QIAGEN, Aarhus, Denmark), assembled *de novo* and compared using MUSCLE⁴⁷. In
470 addition, as a 100-bp gap remained in Peach genome v2.0.a1 in the region immediately upstream of the
471 most likely CG, 25 primers were developed (Table S3) and used for Sanger sequencing of the gap region
472 and the target CG (Genewiz, South Plainfield, NJ, USA). Assembled sequences were then compared
473 and differences between 'S10215' and 'S10216' identified. Sequences of the selected CG and the gap
474 region, were finally analyzed for comparison and possible changes in the coding sequences, as well as
475 changes in the resulting protein, using FGENESH gene-prediction program⁴⁸ with different dicot plant
476 species as model (http://www.softberry.com).
477

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

Gene expression analysis

478 Eight cultivars with contrasted phenotypes and both wild species (Prunus davidiana P1908 and Prunus
479 kansuensis S1429) were selected for expression analysis. Foliar samples were collected from the part of
480 the leaves including leaf-glands, or from the region including the base of the leaf-blade and the upper
481 part of the petiole for the eglandular individuals. Regarding Prunus kansuensis, two samples were
482 collected in order to make comparisons: one from the margin of the leaf-blade where reniform glands
483 were visible, the other from the base, close to the petiole. Samples were immediately frozen in liquid
484 nitrogen. Total RNA was isolated using the Macherey-Nagel® NucleoSpin® RNA Plant kit (Thermo
485 Fischer Scientific, Waltham, MA USA) following manufacturer's instructions. RNA concentration and
486 quality were assessed using a NanoDrop™ ND-1000 spectrophotometer (Thermo Fisher Scientific) and
487 the Agilent 2100 Bioanalyzer System (Agilent, Santa Clara, CA USA). For reverse transcription
488 analysis, primer pairs composed of primers on either side of the second intron of Prupe.7G121100 (Fig.
S3) were designed using Primer3⁴³ and GenScript® Tool (Table S5). One microgram of total RNA per sample was then subjected to cDNA synthesis using the AffinityScript RT kit (Agilent) according to the manufacturer's instructions. A SYBR green real-time PCR assay was thereby carried out in a final volume of 15 µl of a reaction mixture containing 7 µl of 2x Brilliant III SYBR® Green qPCR Master mix (Agilent), 0.5 µM of each primer and 100 ng of cDNA template. Reaction mixtures without cDNA were used as negative controls. Amplification reactions were run in a 96 well plate on a Stratagene Mx3005P (Agilent) under the following conditions: 95°C for 30 s, followed by 40 cycles of denaturation at 95°C for 10 s, annealing at 60°C for 30 s and extension at 72°C for 15 s. Reactions were performed using four biological and three technical replicates for each sample. Amplification values were then normalized using two genes as constitutive controls, as recommended by Bustin et al.⁴⁸: *PpTEF2* (translation elongation factor 2) and *PpRPL13* (60S ribosomal protein L13), both having previously been tested and selected for their stability. Two-way analysis of variance (ANOVA) was used to assess the independent effect of the presence of EFNs and that of the insertion on the expression of *Prupe.7G121100*. Tests were performed using a script of 'RqPCRAnalysis' R-package⁴⁹ customized to generate box-plots with R studio⁵⁰. Significance threshold was set to *p* <0.01.

## 3' RACE PCR

Transcripts of reniform and eglandular accessions were amplified using the Invitrogen 3' Rapid Amplification of cDNA Ends (3' RACE) system (Thermo Fisher Scientific). The 3' RACE procedure was carried out as recommended by the supplier. The first strand cDNA was synthesized using 1 µg of total RNA of each of the individuals and the adapter primer (AP) targeting the poly(A) region of the mRNA. The synthesis reaction was followed by the amplification of the target cDNA in a final volume of 50 µl containing 2 µl (1/10) of the above reaction, 1x reaction buffer, 0.2 mM each dNTP, 1.5 mM MgCl2, 0.2 µM each of the following primers, the antisense abridged universal amplification primer (AUAP) provided in the kit and a custom sense primer developed in the second exon of the gene (Table S6), and 2.5 U of GoTaq® Hot Start Polymerase (Promega). Amplification reactions were run on an Eppendorf Mastercycler epgradient (Eppendorf AG, Hamburg, Germany) under the following conditions: 94°C (2min) followed by 35 cycles at 94°C (45 sec), 57°C (45 sec), 72°C (1.5 min) and a
final extension at 72°C (5 min). PCR products were visualized in a 1.5 % agarose gel stained with ethidium bromide. Nested amplification reactions were then carried out as above using 1/5 of the second reaction, the antisense AUAP primer and additional custom sense primers developed downstream of the first sense primer. In addition, as numerous stretches of poly(A) were included in the sequence of the transposon, which interfered in the hybridization of the adapter primer (AP) to the poly(A) region of the mRNA, antisense primers were developed based on each of the predictions derived from the analysis with FGENESH (Table S6). PCRs were carried out as above except for the annealing temperature which was lowered to 55°C. Resulting amplicons were then sent for sequencing (Genewiz). Finally, upstream regions were amplified and sequenced to obtain the complete transcripts.

## Development of the diagnostic marker ASPP900

One primer-triplet based on the PCR-based KASP™ method (https://www.biosearchtech.com/) was developed in order to differentiate each of the three phenotypes encountered (Table S4). It was composed of one forward primer specific of the glandular phenotype (20 nucleotides in CDS 3 of *Prupe.7G121100* starting 15 nucleotides before the insertion position), one forward primer specific of the eglandular phenotype (18 nucleotides astride the 9 last nucleotides of the transposon and the first 9 nucleotides of CDS3 after the insertion), and one common reverse primer (20 nucleotides in CDS3, starting 75 nucleotides downstream of the insertion point). Positions of the primers on the sequence are shown in Fig. S1.

## Acknowledgements

We thank Amélie Emanuel currently at the UMR-BPMP (Montpellier) and Céline Roques for their involvement in the initial genotyping and phenotyping work of the mapping population. We are grateful to Henri Duval from GAFL, for having managed the PeachReseq project with MGX. We thank Christophe Tuero from GAFL, for his help in the phenotyping and Marine Delmas from the *Prunus-Juglans* Genetic Resource Centre, Domaine des Jarres, 33210 Toulenne (France) for providing peach samples. We also thank Caroline Le Baron for her help with 3' RACE PCR and we are grateful to the technical staff of the experimental domains of 'Saint Maurice' and 'Les Pins de l'Amarine' (INRAE-
542 UGAFL) for their technical contribution to tree management of the peach collection and the population

543 5392². Finally we thank Jean-Luc Gallois (UGAFL) for his help and suggestions as well as Jean-Luc

544 Poëssel (UGAFL) for proofreading the manuscript.

545 **Conflict of interests**

546 The authors declare no potential conflict of interests of any kind.

547 **Data Availability**

548 The datasets supporting the current study are available from the corresponding author upon request, in

549 strict accordance with the policy of the Journal.

550 **Supplementary information** accompanies the manuscript on the *Horticulture Research* website

551 http://www.nature.com/hortres

552 **References**

553 1 Gregory, C.T. The taxonomic value and structure of the peach leaf glands. *N.Y. Cornell Agric. Exp.*

555 *Sta. Bull. 365*, 183-222 (1915).

556 2 Okie, W.R. Handbook of peach and nectarine varieties. USDA-ARS, Byron, GA (1998).

557 3 Koptur, S. Extrafloral nectary-mediated interactions between insects and plants. Insect-Plant

558 Interactions. Vol. IV (ed. E.A. Bernays), pp. 81-129. CRC Press, Boca Raton, FL, USA (1992).

559 4 Marazzi, B., Bronstein, J.L. & Koptur, S. The diversity, ecology and evolution of extrafloral nectaries:

560 current perspectives and future challenges. *Ann. Bot.* **111**, 1243-1250 (2013).

561 5 Mathews, C.R., Bottrell, D.G. & Brown, M.W. Extrafloral nectaries alter arthropod community

562 structure and mediate peach (Prunus persica) plant defense. *Ecological Applications* **19**, 722-

563 730 (2009).

564 6 Jones, I.M., Koptur, S. & von Wettberg, E.J. The use of extrafloral nectar in pest management:

565 overcoming context dependence. *J. of Appl. Ecol.* **54**, 489-499 (2017).
566 7 Weber, M.G., Porturas, L.D. & Taylor, S.A. Foliar nectar enhances plant-mite mutualisms: the effect of leaf sugar on the control of powdery mildew by domatia-inhabiting mites. *Ann. Bot.* **118**(3), 459-466 (2016).

569 8 Scorza, R. & Sherman, W.B. Peaches. pp 325-440 In: Janick J, Moore JN (eds.) Fruit breeding, Volume 1: Tree and tropical fruits. John Wiley & Sons, NY (1996).

571 9 Watking, W. & Brown A.G. Genetic response to selection in cultivated plants: Gene frequencies in varieties of *Prunus persica*. *Proc. Soc. Lond. B. Biol Sci.* **145**, 337-347 (1956).

573 10 Saunier, R. Contribution to the study of relationships between certain characteristics of simple genetic determination in the peach tree and susceptibility of peach cultivars to oïdium, *Sphaerotheca pannosa* (Wallr.). *Lev. Ann. Amelior. Plant* **441**, 235-243 (1973).

576 11 Weinhold, A.R. The orchard development of peach powdery mildew. *Phytopathology* **51**, 478-481 (1961).

578 12 Pascal, T., Pfeiffer, F. & Kervella, J. Powdery mildew in the peach cultivar Pamirkij 5 is genetically linked to the *Gr* gene for leaf color. *Hort. Sci.* **45**(1), 150-152 (2010).

580 13 Connors, C.H. Inheritance of foliar glands of the peach. *Proc. Am. Soc. Hortic. Sci.* **18**, 21-27 (1921).

581 14 Dettori M.T., Quarta, R. & Verde, I. A peach linkage map integrating RFLPs, SSRs, RAPDs, and morphological markers. *Genome* **44**, 783-90 (2001).

583 15 Dirlewanger, E. et al. Comparative mapping and marker-assisted selection in Rosaceae fruit crops. *Proc. Natl. Acad. Sci. USA* **101**, 9891-9896 (2004).

585 16 Verde, I., Quarta, R., Cedrola, C. & Dettori, M.T. QTL analysis of agronomic traits in a BC1 peach population. *Acta Hortic.* **592**, 291-297 (2002).

587 17 Dabov, S. Malo Konare - a new canning peach variety. *Rastenievj Nauki* (1985).
588

18 Micheletti, D. et al. Whole-genome analysis of diversity and SNP-major gene association in peach germplasm. *PloS one* **10**(9), e0136803. https://doi.org/10.1371/journal.pone.0136803

589

590

591

19 Verde, I. et al. Development and evaluation of a 9 K Array for peach by internationally coordinated SNP detection and validation in breeding germplasm. *PloS one* **7**(4), e35668. https://doi.org/10.1371/journal.pone.0035668

592

593

594

20 Wang, L. et al. Evolutionary origin of Rosaceae-specific active non-autonomous hAT elements and their contribution to gene regulation and genomic structural variation. *Plant Mol. Biol.* **91** (1-2), 179-191 (2016).

595

596

597

21 Bennetzen, J.L. & Wang, H. The contributions of transposable elements to the structure, function, and evolution of plant genomes. *Annu. Rev. Plant Biol.* **65**, 505-530 (2014).

598

599

22 Chen, J., Hu Q., Zhang, Y., Lu C & Kuang, H. P-Mite : a database for plant miniature inverted-repeat transposable elements. *Nucleic Acids Res.* **42**, Database issue (2013).

600

601

23 Lu, C., Chen J., Zhang, Y., Hu, Q., Su W & Kuang, H. Miniature inverted-repeat transposable elements (MITEs) have been accumulated through amplification bursts and play important roles in gene expression and species diversity in *Oryza sativa*. *Mol. Biol. Evol.* **29**, 1005-1017 (2012).

602

603

604

605

24 Henriksson, E. et al. Homeodomain leucine zipper class I genes in Arabidopsis. Expression patterns and phylogenetic relationships. *Plant Physiology* **139**, 509-518 (2005).

606

607

25 Saddic, L.A. et al. The LEAFY target LMI1 is a meristem identity regulator and acts together with LEAFY to regulate expression of *CAULIFLOWER*. *Development* **133**, 1673-1682 (2006).

608

609

26 Mangan, S. & Alon, U. Structure and function of the feed-forward loop network motif. *Proc. Natl. Acad. Sci USA* **100**(21), 11980-11985 (2003).

610
611 27 Xu, M. et al. Arabidopsis BLADE-ON-PETIOLE1 and 2 promote floral meristem fate and determinacy in a previously undefined pathway targeting APETALA1 and AGAMOUS-LIKE. *The Plant Journal* **63**, 974–989 (2010).

614 28 Vlad, D. et al. Leaf shape evolution through duplication, regulatory diversification, and loss of a homeobox gene. *Science* **343**, 780–787 (2014).

616 29 Vuolo, F. et al. LMI1 homeodomain protein regulates organ proportions by spatial modulation of endoreduplication. *Genes Dev.* **33**, 377 (2018).

618 30 Andres, R.J. et al. Modifications to a LATE MERISTEM IDENTITY1 gene are responsible for the major leaf shapes of Upland cotton (*Gossypium hirsutum* L.). *Proc. Natl. Acad. Sci USA* **114**(1), E57-E66 (2016).

621 31 Maugarny-Calès, A. & Laufs, P. Getting leaves into shape: a molecular, cellular, environmental and evolutionary view. *Development* 145, dev161646 (2018).

623 32 Chang, L.J., Mei, G.F., Hu, Y., Deng, J.Q. & Zhang, T.Z. LMI1-like and KNOX1 genes coordinately regulate plant leaf development in dicotyledons. *Plant Mol. Biol.* **99**, 449-460 (2019).

625 33 Sicard, A. et al. Repeated evolutionary changes of leaf morphology caused by mutations of a homeobox gene. *Curr. Biol.* **24**(16), 1880-1886 (2014).

627 34 Zhao, J.L. et al. Transcriptome analysis in *Cucumis sativus* identifies genes involved in multicellular trichome development. *Genomics* **105**, 296–303 (2015).

629 35 Hu, W. et al. Genetic and evolution analysis of extrafloral nectary in cotton. *Plant Biotech. Jour.* **18** (10), 2081-2095 (2020).

631 36 Bowman, J.L. & Smyth, D.R. CRABS CLAW, a gene that regulates carpel and nectary development in *Arabidopsis*, encodes a novel protein with zinc finger and helix-loop-helix domains. *Development* **126**, 2387–2396 (1999).
634

37 Gross, T., Broholm, S. & Becker, A. CRABS CLAW Acts as a Bifunctional Transcription Factor in Flower Development. *Front Plant Sci.* **9**, 835 (2018).

635

636

38 Lee, J.Y. et al. Activation of *CRABS CLAW* in the nectaries and carpels of Arabidopsis. *Plant Cell* **17**, 25–36 (2005a).

637

638

39 Lee, J.Y. et al. Recruitment of CRABS CLAW to promote nectary development within the eudicot clade. *Development* **132**, 5021–5032 (2005b).

639

640

40 Bhosale, R., Maere, S. & De Veylder, L. Endoreplication as a potential driver of cell wall modifications. *Curr. Opin. in Plant Biol.* **51**, 58-65 (2019).

641

642

41 Desnoues, E. et al. Dynamic QTLs for sugars and enzyme capacities provide an overview of genetic control of sugar metabolism during peach fruit development. *J. Exp. Bot.* **67**(11), 3419–3431 (2016).

643

644

645

42 Thorvaldsdóttir, H., Robinson, J.T. & Mesirov, J.P. Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration. *Briefings in Bioinformatics* **14**(2), 178-192 (2013).

646

647

648

43 Untergasser, A. et al. Primer3-new capabilities and interfaces. *Nucleic Acids Res.* **40**(15), e115 (2012).

649

650

44 Verde, I. et al. The Peach v2.0 release: high-resolution linkage mapping and deep resequencing improve chromosome-scale assembly and contiguity. *BMC Genomics* **18**, 225 (2017).

651

652

45 Van Ooijen, J.W. JoinMap® 4.1, Software for the calculation of genetic linkage maps in experimental populations of diploid species. Kyazma B.V., Wageningen, Netherlands (2012).

653

654

46 Kosambi, D.D. The estimation of map distance from recombination values. *Ann. Eugenics* **12**, 172-175 (1944).

655

656

47 Madeira, F. et al. The EMBL-EBI search and sequence analysis tools APIs in 2019. *Nucleic Acids Res.* **47**(W1), 636–641 (2019).

657
658

48 Solovyev, V., Kosarev, P., Seledsov, I., Vorobyev, D. Automatic annotation of eukaryotic genes, pseudogenes and promoters. *Genome Biol.* 7, Suppl 1: P. 10.1-10.12 (2006).

659

660

49 Bustin, S.A. et al. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. *Clin. Chem.* **55**, 611-622 (2009).

661

662

50 Hilliou, F. & Tran, T. RqPCRAnalysis: Analysis of Quantitative Real-time PCR Data. *Proceedings of the International Conference on Bioinformatics Models, Methods and Algorithms* pages 202-211 (2013).

663

664

665

51 R Core Team: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. (2020).

666