Julca et al. BMC Biology (2020) 18:148 https://doi.org/10.1186/s12915-020-00881-6 RESEARCH ARTICLE Open Access Genomic evidence for recurrent genetic admixture during the domestication of Mediterranean olive trees (Olea europaea L.) Irene Julca1,2,3, Marina Marcet-Houben1,2,4, Fernando Cruz5, Jèssica Gómez-Garrido5, Brandon S. Gaut6, Concepción M. Díez7, Ivo G. Gut2,5, Tyler S. Alioto2,5, Pablo Vargas8 and Toni Gabaldón1,2,4,9* Abstract Background: Olive tree (Olea europaea L. subsp. europaea, Oleaceae) has been the most emblematic perennial crop for Mediterranean countries since its domestication around 6000 years ago in the Levant. Two taxonomic varieties are currently recognized: cultivated (var. europaea) and wild (var. sylvestris) trees. However, it remains unclear whether olive cultivars derive from a single initial domestication event followed by secondary diversification, or whether cultivated lineages are the result of more than a single, independent primary domestication event. To shed light into the recent evolution and domestication of the olive tree, here we analyze a group of newly sequenced and available genomes using a phylogenomics and population genomics framework. Results: We improved the assembly and annotation of the reference genome, newly sequenced the genomes of twelve individuals: ten var. europaea, one var. sylvestris, and one outgroup taxon (subsp. cuspidata)—and assembled a dataset comprising whole genome data from 46 var. europaea and 10 var. sylvestris. Phylogenomic and population structure analyses support a continuous process of olive tree domestication, involving a major domestication event, followed by recurrent independent genetic admixture events with wild populations across the Mediterranean Basin. Cultivated olives exhibit only slightly lower levels of genetic diversity than wild forms, which can be partially explained by the occurrence of a mild population bottleneck 3000–14,000 years ago during the primary domestication period, followed by recurrent introgression from wild populations. Genes associated with stress response and developmental processes were positively selected in cultivars, but we did not find evidence that genes involved in fruit size or oil content were under positive selection. This suggests that complex selective processes other than directional selection of a few genes are in place. Conclusions: Altogether, our results suggest that a primary domestication area in the eastern Mediterranean basin was followed by numerous secondary events across most countries of southern Europe and northern Africa, often involving genetic admixture with genetically rich wild populations, particularly from the western Mediterranean Basin. Keywords: Admixture, Domestication, Genome, Introgression, Olive * Correspondence: toni.gabaldon.bcn@gmail.com 1 Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology, Dr. Aiguader 88, 08003 Barcelona, Spain 2 Department of Experimental and Health Sciences, Universitat Pompeu Fabra (UPF), 08003 Barcelona, Spain Full list of author information is available at the end of the article © The Author(s). 2020 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Julca et al. BMC Biology (2020) 18:148 Page 2 of 25 Background group (southern Spain and Portugal) retained the finger- The Mediterranean olive tree (Olea europaea L. subsp. print of a genetic bottleneck and, surprisingly, was europaea, Oleaceae) is one of the earliest cultivated fruit closely related to cultivated accessions from the Levant trees of the Mediterranean Basin (MB). Current classifi- [12]. Second, the central MB group, which also included cations recognize two taxonomic varieties of O. euro- the cultivars from eastern Spain (Catalonia, Valencia, paea: var. sylvestris (hereafter sylvestris, also named and Balearic Islands), showed signals of recent and ex- oleaster) for wild populations and var. europaea (here- tensive admixture with local wild populations and rela- after europaea) for cultivated forms [1, 2]. Both varieties tively high plastid diversity compared to the other are predominantly out-crossing and have long lifespans, groups of cultivars. These differential patterns, along including a long juvenile phase that can last up to 15 with the fact that many cultivars from the central MB years in natural conditions. The natural distribution of retain wild-like phenotypic characteristics, opened the the Mediterranean olive encompasses all countries of controversial question of a possible minor domestication the MB, although a few wild populations have also been center for olives in the central MB [10, 14]. Recently, the found in northern areas with low occurrence of frost [3]. analysis of the genome of a set of traditional olive culti- Cultivars have historically been planted nearby wild pop- vars and carbonized pit remains from Spanish Roman ulations since ancient times, where they exchange pollen archeological sites and has favored this latter hypothesis. that have resulted in effective crop production and his- However, the phylogenomic and population genomic torical hybridization [4]. analyses applied ignored plastid lineages and some other There is a large body of evidence pointing to the east- important features of domestication [16]. ern MB as the cradle of the first domestication event of Approximate Bayesian computing (ABC) models have olives. According to archeological, palaeobotanical, and been applied to infer the demographic history of olives genetic studies, the crop was domesticated from eastern and were consistent with a primary domestication event wild progenitors around 6000 years ago [5–8]. Once su- in the East [8, 12]. These models also highlighted the perior individuals were selected, clonal propagation paramount role of admixture to account for the diversity made their multiplication and the spread of valuable of the crop. This feature was particularly predominant in agronomic traits possible. Clonal propagation also facili- cultivars from the central MB, where ~ 20% of the gen- tated the spread of cultivars from the eastern to the etic diversity of olives may have been acquired via intro- western MB via classical civilizations such as Phoeni- gression with local wild populations [12]. So far, genetic cians, Greeks, or Romans. After six millennia, olive do- and archeological sources of evidence have agreed with mestication has resulted in a vast number of cultivars of the existence of a major center of domestication for ol- uncertain pedigree that are often geographically re- ives in the Levant, but questions remain about the extent stricted [9]. of admixture and the potential for secondary centers of It remains unclear whether olive cultivars derived from domestication elsewhere in the MB. To gain novel in- a single initial domestication event in the Levant, sights into this open question and into the most recent followed by secondary diversification [8, 10, 11], or evolution of the cultivated olive tree, we sequenced whether cultivated lineages are the result of more than a twelve accessions, including ten representative cultivars, single, independent primary domestication event [7, 12– one wild individual of var. sylvestris, and one individual 16]. Previous studies based on plastid and nuclear from O. europaea subsp. cuspidata to be used as a dis- markers have suggested controversial but not necessarily tant outgroup in our analysis (Table 1). The ten cultivars incompatible domestication scenarios. The reconstruc- were carefully selected to capture the signal of the diver- tion of plastid lineages has yielded unresolved phyloge- gent genetic clusters reported by previous studies among nies due to low plastid diversity [8, 17, 18]. For instance, the cultivated olive germplasm [8, 12]. Together, they 90% of olive cultivars across the MB share the same represent the following: (i) the three genetic pools iden- “eastern like” plastid haplotype [18]. This general result tified in the MB (West: “Picual” and “Lechin de Gran- is congruent with archeological data [5], which suggests ada;” Central: “Arbequina,” “Frantoio,” “Koroneiki,” and a major domestication event in the Levant, possibly “Megaritiki;” East: “Beladi” and “Sorani”); (ii) the main followed by recurrent admixture events with local wild plastid lineages found in the cultivated olive (“Lechin de olives that would have contributed to the crop diversity. Sevilla” and “Megaritiki”—E2.3 and E2.2, respectively; However, nuclear markers showed a more complex “Chemlal de Kabilye”—E3.2; the rest of cultivars—E1.1); pattern. Olive cultivars clustered together into three dif- and (iii) the main cultivar diversity from the most im- ferent gene pools, with a rough geographical corres- portant areas of Mediterranean olive cultivation. For in- pondence to the eastern, central, and western MB [12, stance, “Arbequina” is the most international cultivar 19–22]. The relationship between these groups also due to its adaptation to high-density planting designs, showed two interesting features [12]. First, the western “Picual” covers more than 1.5 million of hectares in Julca et al. BMC Biology (2020) 18:148 Page 3 of 25 southern Spain, “Frantoio” and “Koroneiki” are the pri- than the recently published assembly for cv. Picual (410, mary cultivars in Italy and Greece, respectively, and 451) [16], but lower than the assembly for O. europaea “Beladi” and “Sorani” are widely used in the Levant. Of var. sylvestris (12,567,911) [24]. Nevertheless, Oe9 is note, our selection includes only authenticated cultivars much less fragmented than the sylvestris genome, dis- preserved in the germplasm collection of Cordoba Uni- playing a ~ 5-fold reduction in the number of scaffolds versity—i.e., those that matched by DNA and morpho- (9751 vs 41,261 scaffolds), and has larger scaffolds than logical markers with authentic control samples from its the cv. Picual genome (Additional file 1: Fig. S1 and natural area of origin and/or cultivation [9]. For sylves- Additional file 2: Table S1). The comparison between tris, previous phylogenetic results and field experience Oe9 and sylvestris assemblies (Additional file 1: Fig. S1b) led us to choose a wild individual from an isolated area shows a high level of conserved synteny and many re- in order to avoid potential feral or highly introgressed gions appeared duplicated between the two genomes. trees. Hence, we sampled a tree from a location by the These duplicated syntenic regions showed different coast near the Cantabrian mountains (northern Spain) levels of divergence, with Ks values peaking around me- where the olive tree has not been historically cultivated dian values of 0.32 and 59.19. This supports the occur- and that is 200 km distant from current plantations. rence of more than one ancient polyploidization event in These populations had previously been screened using O. europaea, as proposed earlier [26]. fingerprinting techniques [3] and Sanger sequencing Additionally, we improved the genome annotation in [18]. Oe9 genome, by extended automated functional annota- These twelve newly obtained sequences plus 43 add- tion and by manual curation of some of the genes. Con- itional sequences published in a recent study [16] com- served gene completeness results for both Oe6 and Oe9 plement the three available genomes for the species: the support that the decontamination step (see “Methods”) cv. Farga from eastern Spain, for which we provide here did not result in the loss of sequences coding for import- an improved assembly by anchoring it to a genetic map; ant plant genes. Moreover, the Oe9 genome assembly the cv. Picual; and an oleaster (var. sylvestris) from has clearly higher gene completeness than the sylvestris Turkey [16, 23, 24]. These datasets combined constitute reference (93% vs 87%). The cvs. Picual and Farga (Oe9) the largest and most comprehensive ensemble of olive assemblies have very similar gene completeness values, whole genomic information, which we analyzed here although the percentage of duplicated sequence is twice under a phylogenomics and population genomics frame- in cv. Picual (Additional file 2: Table S1). To better work to shed light on the recent evolution and domesti- characterize duplicated sequence differences between cation of the olive tree. In particular, we addressed the the cvs. Farga (Oe9) and Picual references, two RNAseq question whether genome-wide variation data can disen- libraries (SRR6003535, ERR1406351) were aligned to the tangle scenarios of one versus multiple centers of do- two genome assemblies. Results confirmed that both as- mestication. Additionally, we were interested in finding semblies are very similar in terms of completeness but out whether genetic introgression between wild and cul- show a high percentage of multi-mappings in the cv. tivated trees have historically played a role in the domes- Picual reference (Additional file 2: Table S2), which sug- tication process, to assess the presence of potential gests the presence of some artefactual duplications in domestication bottlenecks and to identify genes and gen- the Picual reference. omic regions under selection. Finally, sequencing of nu- The Oe9 assembly has 4911 more annotated genes clear genomes allowed testing earlier suggestions of than the var. sylvestris genome. A comparison of gene close relationship between cultivars of distant locations sets using BLASTN [27] with identity > 80% and e-value such as southern Spain and the Levant [12]. < 1e−5 cutoffs shows that 5245 Oe9 annotated genes do not have a match in sylvestris, and conversely, 2620 Results genes of sylvestris do not have a match in europaea New assembly version of the reference olive (cultivar (Oe9). Distinct genome annotation methods could partly Farga) genome explain these differences [23, 24]. To have an We improved the available O. europaea var. europaea annotation-independent measure of the differences be- (cv. Farga) genome assembly (Oe6 version) [23] by an- tween both assemblies, we mapped raw reads of both choring it to chromosomes using a publicly available genome projects to the alternative assembly and assessed genetic map [25] and removing 201 scaffolds which coverage of the putative unique genes (see “Methods”). likely represent contaminating sequences (see We first filtered out the genes that have at least 50% of “Methods”). In the final assembly (Oe9), 520.5 Mbp their length with a read coverage higher than 20, which (39.5%) of Oe6 sequence was anchored to 23 linkage resulted in 2115 and 280 unique genes for europaea and groups, of which 288 Mbp (21.8%) were oriented. This sylvestris, respectively. Even when lowering the coverage anchored assembly (Oe9) has a higher N50 (734,380) threshold to 5, europaea still had more unique genes Table 1 O. europaea genomes used in the analysis. The columns show the sample origin, analysis in which it was used, plastid group [17], nuclear group [12], ploidy level, and the source of the data. *Cultivars that are duplicated and are excluded from the main analysis (see “Methods”) O. europaea subsp. Origin Analysis Plastid Nuclear Ploidy Source group group level europaea var. europaea cv. “Farga” Spain (Boadilla/La Senia) Nuclear, plastid, mitochondrial E3.1 Central MB (Q2) 2x ENA (PRJEB4992) europaea var. europaea cv. “Arbequina” Spain (BGMO_UCO) Nuclear, plastid, mitochondrial E1.1 Central MB (Q2) 2x This study Julca et al. BMC Biology europaea var. europaea cv. “Picual” Spain (BGMO_UCO) Nuclear, plastid, mitochondrial E1.1 Western MB (Q1) 2x This study europaea var. europaea cv. “Beladi” Lebanon (BGMO_UCO) Nuclear, plastid, mitochondrial E1.1 Eastern MB (Q3) 2x This study europaea var. europaea cv. “Sorani” Syria (BGMO_UCO) Nuclear, plastid, mitochondrial E1.1 Eastern MB (Q3) 2x This study europaea var. europaea cv. “Manzanilla” Spain Plastid E1.1 – 2x NCBI (FN996972) (2020) 18:148 europaea var. europaea cv. “Bianchera” Italy Plastid E1.1 – 2x NCBI (NC_013707) europaea var. europaea cv. “Chemlal de Kabilye” Algeria (BGMO_UCO) Nuclear, plastid, mitochondrial E3.2 Mosaic (Q2 + Q1) 2x This study europaea var. europaea cv. “Megaritiki” Greece (BGMO_UCO) Nuclear, plastid, mitochondrial E2.2 Central (Q2) 2x This study europaea var. europaea cv. “Lechin de Sevilla” Spain (BGMO_UCO) Nuclear, plastid, mitochondrial E2.3 Mosaic (Q1 + Q2) 2x This study europaea var. europaea cv. “Lechin de Granada” Spain (BGMO_UCO) Nuclear, plastid, mitochondrial E1.1 Western MB (Q1) 2x This study europaea var. europaea cv. “Frantoio” Italy (BGMO_UCO) Nuclear, plastid, mitochondrial E1.1 Central (Q2) 2x This study europaea var. europaea cv. “Koroneiki” Greece (BGMO_UCO) Nuclear, plastid, mitochondrial E1.1 Central (Q2) 2x This study europaea var. sylvestris “sylvestris-T” Turkey Nuclear, plastid, mitochondrial E1.1 – 2x NCBI (PRJNA417827) europaea var. sylvestris “sylvestris-S” Spain (Pechón) Nuclear, plastid, mitochondrial E3 – 2x This study europaea var. sylvestris “W7R224” Morocco Nuclear, plastid, mitochondrial E3 – 2x SRR9860491 europaea var. sylvestris “W11R37” Spain Nuclear, plastid, mitochondrial E3 – 2x SRR9860496 europaea var. sylvestris “W9R302” Albania Nuclear, plastid, mitochondrial E2 – 2x SRR9860497 europaea var. sylvestris “W8R225” Morocco Nuclear, plastid, mitochondrial E1 – 2x SRR9860498 europaea var. sylvestris “W4R183” Spain Nuclear, plastid, mitochondrial E1 – 2x SRR9860503 europaea var. sylvestris “W3R78” Spain (Menorca island) Nuclear, plastid, mitochondrial E2 – 2x SRR9860505 europaea var. sylvestris “W2R74” Spain Nuclear, plastid, mitochondrial E3 – 2x SRR9860532 europaea var. sylvestris “W1R198” Croatia Nuclear, plastid, mitochondrial E2 – 2x SRR9860533 europaea var. europaea cv. “Myrtolia” Greece Nuclear, plastid, mitochondrial E1 – 2x SRR9860492 europaea var. europaea cv. “Klon-14-1812” Albania Nuclear, plastid, mitochondrial E1 – 2x SRR9860493 europaea var. europaea cv. “Barri” Syria Nuclear, plastid, mitochondrial E1 – 2x SRR9860494 europaea var. europaea cv. “Grappolo” Italy Nuclear, plastid, mitochondrial E1 – 2x SRR9860495 europaea var. europaea cv. “Mari” Iran Nuclear, plastid, mitochondrial E1 – 2x SRR9860499 europaea var. europaea cv. “Lianolia Kerkyras” Greece Nuclear, plastid, mitochondrial E2 – 2x SRR9860500 europaea var. europaea cv. “Piñonera” Spain Nuclear, plastid, mitochondrial E1 – 2x SRR9860501 europaea var. europaea cv. “Mavreya” Greece Nuclear, plastid, mitochondrial E1 – 2x SRR9860502 Page 4 of 25 Table 1 O. europaea genomes used in the analysis. The columns show the sample origin, analysis in which it was used, plastid group [17], nuclear group [12], ploidy level, and the source of the data. *Cultivars that are duplicated and are excluded from the main analysis (see “Methods”) (Continued) O. europaea subsp. Origin Analysis Plastid Nuclear Ploidy Source group group level europaea var. europaea cv. “Royal” Spain Nuclear, plastid, mitochondrial E1 – 2x SRR9860504 europaea var. europaea cv. “Llumeta” Spain Nuclear, plastid, mitochondrial E3 – 2x SRR9860506 Julca et al. BMC Biology europaea var. europaea cv. “Barnea” Israel Nuclear, plastid, mitochondrial E1 – 2x SRR9860507 europaea var. europaea cv. “Majhol-1013” Syria Nuclear, plastid, mitochondrial E1 – 2x SRR9860508 europaea var. europaea cv. “Mastoidis” Greece Nuclear, plastid, mitochondrial E1 – 2x SRR9860509 europaea var. europaea cv. “Abbadi Abou Gabra-842” Syria Nuclear, plastid, mitochondrial E1 – 2x SRR9860510 europaea var. europaea cv. “Fishomi” Iran Nuclear, plastid, mitochondrial E1 – 2x SRR9860511 (2020) 18:148 europaea var. europaea cv. “Abou Satl Mohazam” Syria Nuclear, plastid, mitochondrial E1 – 2x SRR9860512 europaea var. europaea cv. “Forastera de Tortosa” Spain Nuclear, plastid, mitochondrial E3 – 2x SRR9860514 europaea var. europaea cv. “Dokkar” Tunisia Nuclear, plastid, mitochondrial E3 – 2x SRR9860516 europaea var. europaea cv. “Jabali” Syria Nuclear, plastid, mitochondrial E1 – 2x SRR9860517 europaea var. europaea cv. “Majhol-152” Syria Nuclear, plastid, mitochondrial E1 – 2x SRR9860518 europaea var. europaea cv. “Picudo” Spain Nuclear, plastid, mitochondrial E1 – 2x SRR9860519 europaea var. europaea cv. “Maarri” Syria Nuclear, plastid, mitochondrial E1 – 2x SRR9860520 europaea var. europaea cv. “Temprano” Spain Nuclear, plastid, mitochondrial E2 – 2x SRR9860521 europaea var. europaea cv. “Menya” Spain Nuclear, plastid, mitochondrial E2 – 2x SRR9860522 europaea var. europaea cv. “Koroneiki”* Greece Nuclear, plastid, mitochondrial E1.1 – 2x SRR9860523 europaea var. europaea cv. “Manzanilla de Sevilla” Spain Nuclear, plastid, mitochondrial E1 – 2x SRR9860524 europaea var. europaea cv. “Uslu” Turkey Nuclear, plastid, mitochondrial E1 – 2x SRR9860525 europaea var. europaea cv. “Lechin de Sevilla”* Spain Nuclear, plastid, mitochondrial E2.3 – 2x SRR9860526 europaea var. europaea cv. “Ocal” Spain Nuclear, plastid, mitochondrial E1 – 2x SRR9860527 europaea var. europaea cv. “Abou Kanani” Syria Nuclear, plastid, mitochondrial E1 – 2x SRR9860528 europaea var. europaea cv. “Verdial de Velez-Malaga-1” Spain Nuclear, plastid, mitochondrial E1 – 2x SRR9860529 europaea var. europaea cv. “Arbequina”* Spain Nuclear, plastid, mitochondrial E1.1 – 2x SRR9860530 europaea var. europaea cv. “Chemlal de Kabylie”* Algeria Nuclear, plastid, mitochondrial E3.2 – 2x SRR9860531 europaea var. europaea cv. “Leccino” Italy Nuclear, plastid, mitochondrial E1 – 2x SRR9860534 europaea var. europaea cv. “Manzanillera de Huercal Overa” Spain Nuclear, plastid, mitochondrial E2 – 2x SRR9860535 europaea var. europaea cv. “Hojiblanca” Spain Nuclear, plastid, mitochondrial E1 – 2x SRR9860536 europaea var. europaea cv. “Kalamon” Greece Nuclear, plastid, mitochondrial E1 – 2x SRR9860537 europaea var. europaea cv. “Picual”* Spain Nuclear, plastid, mitochondrial E1.1 – 2x SRR9860538 europaea var. europaea cv. “Zarza” Spain Nuclear, plastid, mitochondrial E2 – 2x SRR9860539 Page 5 of 25 Table 1 O. europaea genomes used in the analysis. The columns show the sample origin, analysis in which it was used, plastid group [17], nuclear group [12], ploidy level, and the source of the data. *Cultivars that are duplicated and are excluded from the main analysis (see “Methods”) (Continued) O. europaea subsp. Origin Analysis Plastid Nuclear Ploidy Source group group level europaea var. europaea cv. “Morrut” Spain Nuclear, plastid, mitochondrial E1 – 2x SRR9860540 europaea var. europaea cv. “Frantoio”* Italy Nuclear, plastid, mitochondrial E1.1 – 2x SRR9860541 Julca et al. BMC Biology europaea var. sylvestris “Stavrovouni 11” Cyprus Plastid E1.4 – 2x NCBI (HF558645) europaea var. sylvestris “Haut Atlas 1” Morocco (High Atlas) Plastid E2 – 2x NCBI (NC_015401) europaea var. sylvestris “Gue de Constantine 20” Algeria: Gue de Constantine, Algiers Plastid E3 – 2x NCBI (FN997651) europaea var. sylvestris “Oeiras 1” Portugal Plastid E3.1 – 2x NCBI (MG255763) europaea var. sylvestris “Vallee du Fango 5” France Plastid E2.1 – 2x NCBI (MG255762) (2020) 18:148 maroccana “Immouzzer S1” Morocco (High Atlas) Plastid M – 6x NCBI (NC_015623) guanchica “La Gomera 10” Spain Plastid M-g1 – 2x NCBI (MG255764) laperrinei “Adjelella 10” Algeria Plastid E1-l1 – 2x NCBI (MG255765) cuspidata - R Reunion island Nuclear, plastid, mitochondrial A – 2x This study cuspidata “Almihwit 5.1” Yemen Plastid C2 – 2x NCBI (FN996943) cuspidata “Guanghzou 1” China Plastid C1 – 2x NCBI (FN996944) cuspidata “Maui 1” USA (Hawaii-Maui) Plastid A – 2x NCBI (NC_015604) cuspidata “Menagesha Forest 14” Ethiopia Plastid C2 – 2x NCBI (MG255760) Page 6 of 25 Julca et al. BMC Biology (2020) 18:148 Page 7 of 25 (1756) than sylvestris (102). Of these, we discarded 131 (710,737–769,995 bp) [40]. The coding regions in the olive Oe9-specific genes as possible contaminations as their mitochondrion comprise 46 protein-coding genes, 3 ribo- first BLAST hit fell outside plants. Thus, some of the somal RNA genes, and 26 transfer RNA genes (Additional genes uniquely found in the Oe9 assembly may repre- file 1: Fig. S3, Additional file 2: Table S4). sent true differences in terms of gene content. We assessed the presence of these genes in all the other Contrasting genetic diversity patterns in organellar and samples and found that roughly 50% of those genes were nuclear genomes not found in any of them (1092). We removed 104 Oe9 We used this improved reference genome assembly unique genes that appeared clustered into contigs with (Oe9) to call SNPs at the nuclear, plastid, and mitochon- no other genes and had blast hits in other non-plant ge- drial genomes for all subsp. europaea individuals for nomes, potentially indicating remaining contamination. which whole genome sequence information is available. The remaining Oe9 unique genes are often found in For this, we compiled data comprising a total of 56 clusters within larger contigs, indicating they may have unique subsp. europaea individuals, combining samples been lost together. Among the 1023 genes that are ab- sequenced in this project and samples from recent publi- sent from the sylvestris genome and are found in Oe9 cations [16, 23, 24] (see “Methods”, Table 1). This new and at least one of the other cultivars, there are func- dataset included 46 different cultivars and 10 individuals tions associated with stress response, such as HIPPs described as var. sylvestris. Altogether, for wild and culti- (heavy metal-associated isoprenylated plant proteins) vated olives (subsp. europaea), we obtained a total of 24, [28], LEA (late embryogenesis abundant) [29], and sali- 724,756 polymorphic positions uniformly distributed cylic acid-binding [30]. Other genes are associated with along the nuclear genome (Additional file 1: Fig. S4 and growth and development. This is the case of RALF S5), 85 in the plastid genome, and 3979 in the mito- (rapid alkalinization factor), which has been shown to chondrial genome (see Additional file 1: Fig. S2, S3). In arrest root growth and development in tomato and Ara- the plastid genome, a large region (~ 25 kb) was found bidopsis [31], and caffeoyl shikimate esterase (CSE), to be fully conserved and devoid of SNPs in all analyzed which is an enzyme central to the lignin biosynthetic individuals (Additional file 1: Fig. S2). This region in- pathway [32]. It is well known that lignin biosynthesis cludes the largest plastid gene, ycf2, which has also been contributes to plant growth, tissue and organ develop- found to exhibit low rates of nucleotide substitution in ment, and response to a variety of biotic and abiotic other plants [41]. This gene is essential for plant sur- stresses [33]. Some other Oe9 unique genes were associ- vival; however, its exact function is unknown [42, 43]. ated with seed dormancy and sugar signaling, DOG1 This conserved region also comprises other genes, in- [34], and positive regulation of germination, PELPK1 cluding ycf15, rps7, rpl2, ndhB, rRNA, and tRNA. All in- [35] (see Additional file 2: Table S3). Only two sylvestris dividuals presented similar amounts of nuclear unique genes had annotated functions (Additional file 2: polymorphisms relative to the Oe9 reference (Fig. 1a). Table S3). One corresponds to GSH-induced LITAF, Interestingly, the sylvestris from northern Spain (sylves- which negatively regulates hypersensitive cell death in tris-S) has a higher number of homozygous SNPs and a Arabidopsis [36]. The other one corresponds to FAR1 lower number of heterozygous SNPs (Fig. 1a). The other (far-red-impaired response)-related sequence, with roles sylvestris have different patterns, five of them (W2R74, in diverse developmental and physiological processes W9R302, W3R78, W7R224, W11R37) have a high num- [37, 38]. ber of homozygous SNPs, but similar number of hetero- In addition, we assembled the plastid and mitochon- zygous SNPs as some cultivars, while the other four drial genomes of the cv. Farga, which were not provided sylvestris (W8R225, W4R183, W1R198, sylvestris-T) have as separate assemblies in the previous release [23] (see a number of homozygous and heterozygous SNPs in the “Methods”). The final assembly of the plastid genome range of variation found in cultivars (Fig. 1a). comprised 155,658 base pairs (bp) (Additional file 2: Strikingly, the relative number of organellar polymor- Table S4), in agreement with previously reported olive phisms for the different individuals did not follow the plastid sequences, which range from 155,531 to 155,896 patterns observed in the nuclear genome (Fig. 1b). The bp [17, 39, 40]. We annotated 130 genes out of the 130– plastid and mitochondrial genomes of four sylvestris 133 genes reported for other olive plastid genomes (see (W7R224, W2R74, W11R37, and sylvestris-S) and four Additional file 1: Fig. S2, Additional file 2: Table S4) [17, cvs. (Chemlal de Kabilye, Forastera de Tortosa, Llumeta, 39], of which 85 are protein coding genes, 37 are transfer and Dokkar) show a notably lower number of SNPs rela- RNAs, and eight are ribosomal RNAs. The final assem- tive to the cv. Farga reference genome (Fig. 1b). Se- bly of the mitochondrial genome has a size of 755,572 quence variation in the plastid was remarkably low, in bp (Additional file 2: Table S4), which is similar to that agreement with previous studies [17]. Plastid genomes of previously sequenced wild olive mitochondrial genomes can be arranged into three groups based on nucleotide Julca et al. BMC Biology (2020) 18:148 Page 8 of 25 Fig. 1 (See legend on next page.) Julca et al. BMC Biology (2020) 18:148 Page 9 of 25 (See figure on previous page.) Fig. 1 SNP density (SNPs/kb) in sequenced individuals. a Homozygous versus heterozygous SNPs for each accession, relative to the cv. Farga reference. Dot size correlates with the total amount of SNPs. All the cultivars are marked in green and var. sylvestris in blue. b SNP densities for the plastid and mitochondrial genomes. c Plot showing the relative position and identity of plastid SNPs compared to the cv. Farga (reference genome). Bars on the bottom show the main plastid haplotypes of the individuals as described by Besnard et al. [17, 44] polymorphisms, which are congruent with the three introgression using a model-based genetic structure ana- main plastid lineages (E1, E2, and E3) already described lysis (see “Methods”). A phylogeny based on nuclear (Fig. 1c) [17, 44]. More interestingly, eight individuals, polymorphisms (Fig. 2a) places six sylvestris individuals four cultivars and four oleasters, share the same haplo- (sylvestris-S, W11R37, W7R224, W3R78, W2R74, type as cv. Farga (E3), with 3 or less SNPs of difference W9R302) as ancestral lineages of all the other subsp. (Fig. 1c). Hence, the nuclear and organellar genomes europaea individuals and closer to subsp. cuspidata. The presumably reflect different evolutionary histories. Not- other four sylvestris appear to branch within the culti- ably, the organellar genomes of the cvs. Farga, Chemlal vated individuals, with three of them (W1R198, de Kabilye, Forastera de Tortosa, Llumeta, and Dokkar, W4R183, W8R225) apparently diverging after the cv. and the wild individuals W7R224, W2R74, W11R37, and Dokkar and ancestral to the remaining cultivars. The sylvestris-S, all from the western MB, show a very close other oleaster, sylvestris-T, is intermingled within the genetic relationship (Fig. 1b, c). As organelles are mater- main cultivars from Lebanon and Syria. This raised the nally inherited in olive [45], our results suggest that question of the authenticity or pureness of cv. Dokkar these cultivars and wild individuals share a recent com- and sylvestris-T as cultivated or wild individuals, respect- mon maternal ancestor. Similarly, “Megaritiki” (from ively. Genetic structure analyses including all 56 individ- Greece), “Lianolia Kerkyras” (Greece), “Temprano” uals suggested the presence of two distinct ancestral (Spain), “Menya” (Spain), Manzanillera de Huercal Overa genetic pools, which are differentially present among the (Spain), “Zarza” (Spain), and “Lechin de Sevilla” (Spain) different individuals (Fig. 2b, Additional file 2: Table S5). share the same plastid haplotype (E2) as the sylvestris Based on this, we distinguished three groups: one com- W1R198 (from Croatia), W9R302 (Albania), W3R78 posed only by cluster 1 (four wild individuals: sylvestris- (Spain), which can also be found in wild populations but S, W11R37, W7R224, W3R78), another composed by exclusively from central and western Mediterranean cluster 2 (14 individuals: sylvestris-T, “Uslu,” “Barri,” areas [18]. The other cultivars and the three oleasters “Jabali,” “Majhol-152,” “Maarri,” “Sorani,” “Beladi,” “Fish- (W4R183, W8R225, sylvestris-T) share the E1 haplotype. omi,” “Mari,” “Abou Satl Mohazam,” “Majhol-1013,” Furthermore, we found a very close relationship at the “Abbadi Abou Gabra-842,” “Abou Kanani”), and the lar- organellar and nuclear levels between the eastern culti- gest consisting of a mixture of the two clusters in differ- vars from Syria, Iran, Lebanon, and Turkey (e.g., “Sorani, ent proportions (the remaining 38 individuals). The ” “Jabali,” “Beladi,” “Fishomi,” “Mari,” and “Uslu”) and same three genetic clusters were also supported by a cultivars from southern Spain (e.g., “Picual,” “Lechín de principal component analysis (PCA; Additional file 1: Granada,” “Manzanilla de Sevilla,” “Hojiblanca,” and Fig. S6). When we analyzed the oleaster samples, we ob- “Picudo”). This finding provides support for a previously served that sylvestris-S, W11R37, W7R224, and W3R78 proposed hypothesis of a recent bottleneck affecting a are composed only by cluster 1 (Fig. 2b), which may rep- subset of western cultivars [12]. This “local” bottleneck resent the genetic pool of western MB wild populations. affected only western olive cultivars (mostly from south- The oleasters W2R74, W9R302, W1R198, W4R183, and ern Spain and Portugal), and it was probably related to W8R225 are composed largely by cluster 1 mixed with the introduction of olive germplasm into southern Spain different amounts of cluster 2 (Additional file 2: Table during the Muslim period. This period began c. 700 AD, S5), which could be product of gene flow between wild lasted eight centuries, and possibly reshaped the culti- and cultivated populations. The oleaster sylvestris-T is vated olive germplasm of the Iberian Peninsula due to composed only by cluster 2, similar to cultivars from purported migration of cultivars from the Levant and Lebanon, Iran, Turkey, and Syria. Based on this and the the possible substitution of the previously cultivated results of the nuclear tree, we suggest that sylvestris-T local germplasm [16]. represents a feral individual and may have been misiden- tified. Interestingly, all cultivars shared cluster 2, which Population structure and phylogenetic relationships of is pervasive and enriched in cultivars from the eastern the olive tree MB in different proportions, suggesting that this may be We reconstructed the phylogenetic relationships of wild a consistent fingerprint of the primary domestication and cultivated olives using nuclear, plastid, and mito- event that took place in this area [18]. The dominance chondrial SNPs separately. In addition, we tested of this genetic background among cultivars suggests that Julca et al. BMC Biology (2020) 18:148 Page 10 of 25 Fig. 2 Maximum likelihood species tree derived from the nuclear SNPs data. a Nuclear phylogeny. Cultivated olives are shown in green and wild olives in blue. The geographical location of the accession and the plastid haplotype are indicated. Only bootstrap values below 100% are shown. b Bayesian clustering for the nuclear SNP data estimated in Structure v2.3. Structure bar plot shows the genetic clusters differentiated by color. c Heatmap showing the D-statistic and its p value. Red color indicates higher D-statistics, and more saturated colors indicate greater significance most cultivars mainly derive from a common primary To reconstruct the plastid tree, we included the avail- domestication process. However, the presence of add- able plastid genomes of four additional accessions of O. itional gene pools within the cultivars depicts patterns europaea subsp. cuspidata and one individual of other that could have resulted from preferential selection of subspecies that are publicly available (see “Methods”). In genetic variants among standing variation, from separate this phylogeny (Fig. 3a), cuspidata individuals grouped domestication events, or from introgression events with together in congruence with previous results [17, 44]. wild populations. Remarkably, consistent with the plastid polymorphism To further investigate the genetic distance between patterns discussed earlier (Fig. 1b,c), we observed incon- populations, we calculated the fixation index (Fst) of the gruencies between the organellar and nuclear trees oleaster population composed only of cluster 1 (sylves- (Fig. 2a and Fig. 3), which provide support for a phylo- tris-S, W11R37, W7R224, W3R78) and two groups of genetic signal of hybridization [46–48]. In particular, all cultivars. The first group is composed only of individuals cultivars, except the cv. Dokkar (plus the possible feral with cluster 2 (cultivars_set1, eastern cultivars), and the sylvestris-T), are monophyletic in the nuclear tree, second is composed of individuals showing a mixture of whereas the plastid tree shows that cultivars are poly- clusters 1 and 2 (cultivars_set2, western cultivars). This phyletic, falling into three independent clades, each of analysis shows that the divergence between western ole- them associated with different sylvestris samples. More- aster and mixed cultivars (weighted Fst = 0.24) was over, attending to the position of samples from other smaller than the divergence of the western oleaster and subspecies (such as guanchica, marrocana, and laperri- eastern cultivars (weighted Fst = 0.41). nei), which are intermingled within the subspecies Julca et al. BMC Biology (2020) 18:148 Page 11 of 25 europaea, the plastid genome suggests a polyphyly of the similar pattern has been observed in a previous study subspecies europaea. that combined a large sample of cultivated olives and A particularly strong phylogenetic incongruence ob- oleasters, in which most cultivars were assigned to the served among cultivars involves the cultivars with plastid eastern genetic pool, even those with maternal lineages type E3 (“Farga,” “Dokkar,” “Chemlal de Kabilye,” “Llu- that originated from the western Mediterranean basin meta,” and “Forastera de Tortosa”), which cluster to- [11, 18]. Indeed, the largest wild populations of olives gether with the other cultivated olives in the nuclear are found in the western Mediterranean. In sum, the tree, but are sister to sylvestris accessions with the same plastid tree suggests genetic contributions, at least in the plastid type (E3), far from the other cultivars in both maternal lineage, of three different genetic pools that plastid and mitochondrial trees (Fig. 3a,b). These results may be in agreement with multilocal introgression or suggest that the maternal line of cultivars with plastid domestication processes, whereas the nuclear tree sug- type E3 originated from western Mediterranean wild ol- gests a dominant, congruent signal shared by all se- ives (which carry the sample E3 plastid haplotype), while quenced cultivars that is consistent with a unique the paternal line originated mostly from domesticated primary domestication center in the eastern MB. These individuals from the eastern Mediterranean basin. A contrasting patterns reinforce the idea that cultivars are Fig. 3 Maximum likelihood species tree derived from the organellar SNPs data. a Plastid phylogeny. Cultivated olives are shown in green and wild olives in blue. The geographical location of the accession and the plastid haplotype are indicated. Only bootstrap values below 100% are shown. b Mitochondrial phylogeny. The colors and characteristics are the same as in a Julca et al. BMC Biology (2020) 18:148 Page 12 of 25 either from the eastern genetic pool or admixed forms genomic regions with a D-statistic higher than 0.5 (Add- [6, 10, 18] and support secondary admixture processes itional file 1: Fig. S7a,b, Additional file 2: Table S6). Indi- in the western Mediterranean basin, with contribution viduals with higher D-statistics have larger genomic from western populations of var. sylvestris, as clearly regions with signatures of introgression. Then, we shown by the plastid lineages. assessed whether specific functions were enriched in genes located within introgressed regions or non- Genetic introgression from western var. sylvestris genetic introgressed regions (common to all five cultivars) and pools among cultivated olives found no significant enrichment in either case. Previous phylogenetic analyses have suggested the exist- In agreement with the split tree, the individuals with ence of ancestral inter-subspecies hybridization in deep the genetic cluster 2 (Fig. 2b) show no signs of introgres- nodes of the evolutionary tree of O. europaea [26]. Our sion with any of the four sylvestris with genetic cluster 1. population genomics results discussed above and a re- This contrasts with the other analyzed cultivars that cently published transcriptome-based population ana- have been strongly introgressed from wild olives of the lysis [11] suggest recent intraspecific genetic admixture western MB (here represented by sylvestris-S, W7R224, between western cultivars and western sylvestris. To in- W11R37, and W3R78). Importantly, the level of intro- corporate genetic admixture into a phylogenetic frame- gression is largely independent of the plastid haplotype work, we reconstructed a split network tree based on (all cultivars of plastid haplotypes E2 and E3 and some nuclear genome data (Fig. 4), which revealed a heavily of E1 show signs of introgression), suggesting that, in at reticulated structure mostly affecting the relationships least some cases, different introgressions have occurred among all olive samples. In particular, most samples in the maternal and paternal lineages, independently. with the exception of oleasters sylvestris-S, W11R37, Most of the introgressed cultivars originate from areas W7R224, and W3R78 appear in a heavily reticulated of the western or central MB. The 14 non-introgressed area. Consistent with the population genomics results, genotypes were all sampled from the eastern MB sylvestris-T appears well embedded within cultivars and (Lebanon, Syria, Iran, and Turkey), close to the pur- together with other 13 cultivars show a distant connec- ported origin of olive primary domestication. Overall, tion from the other accessions. Also, cv. Dokkar and our genome-based results are consistent with earlier re- three other oleasters appear in the reticulated portion of sults on a broader dataset and on transcriptome-based the tree, but further apart from the other individuals. analysis [11]. Overall, such a reticulated network can be the result of hybridization (including introgression) or incomplete Demographic analysis suggests a population bottleneck lineage sorting. In order to detect evidence of introgres- coupled to the early domestication period sion among the olive samples, we ran the ABBA-BABA To investigate the differences in genetic diversity be- test [49] using the program Dsuite v0.1 r3 [50] (see tween cultivars and wild individuals, we used the three “Methods”) for all quartets of the phylogenetic tree defined populations. The first group consisted of four (Fig. 2a). This analysis supports the hypothesis that oleasters with only genetic cluster 1, the second of culti- introgression occurred between the oleasters with a vars with genetic cluster 2, and the third of cultivars unique genetic background (Fig. 2b: sylvestris-S, with mixed genetic clusters. For the last group, we ex- W7R224, W11R37, and W3R78) and eight individuals cluded cv. Dokkar for its phylogenetic position (Fig. 2b) (five cultivars: “Dokkar,” “Chemlal de Kabilye,” “Megari- (see “Methods”). The average pairwise nucleotide diver- tiki,” “Lianolia Kerkyras,” “Menya,” and three oleaster: sity based on 20-kb windows in the population of sylves- W4R183, W8R225, W1R198) with a significant D-statis- tris (3.66 × 10− 3) is slightly higher than that in cultivars_ tic > 0.43 (Fig. 2c). Interestingly, the cv. Dokkar shows set1 (3.52 × 10− 3) and lower than that in cultivars_set2 the highest level of introgression (D-statistic > 0.53), (4.26 × 10− 3). A lower nucleotide diversity in cultivars followed by the three wild individuals. These oleasters when compared with oleaster has been observed in pre- (W4R183, W8R225, W1R198) also appear closer to vious studies based on inter-simple sequence repeats other cultivars in the nuclear tree and the split network (ISSRs) [3], allozyme polymorphisms [51], simple se- tree (Fig. 2a and Fig. 4), which may point to a possible quence repeats (SSRs) [12], and plastid DNA variation feral origin. To study whether the D-statistic is homoge- [17]. Similarly, a recent transcriptome-based analysis re- neous or variable along the genome, we investigated the ported slightly lower genetic diversity in cultivars com- D-statistic over sliding windows (see “Methods”) using pared to wild olives, leading to the suggestion of a weak- the five cultivars with higher signals of introgression moderate population bottleneck during domestication (“Dokkar,” “Chemlal de Kabilye,” “Megaritiki,” “Lianolia [11]. Kerkyras,” “Menya”). This analysis shows that the five in- In general, lower genetic diversity in cultivars is com- dividuals have different amounts of introgressed monly associated with genetic bottlenecks during Julca et al. BMC Biology (2020) 18:148 Page 13 of 25 Fig. 4 SplitsTree derived from nuclear SNPs. All the cultivars are marked in green and var. sylvestris in blue. The neighbor net method is used here to explore data conflict and not to estimate phylogeny domestication [52]. The difference between wild and tool that handles unphased genomes [55] (see cultivated trees observed here is less pronounced than “Methods”). The results of this analysis show evidence that of many domesticated plant species [52, 53]. How- for a continuous decline in population size starting ap- ever, perennial crops often do not show evident domesti- proximately ~ 14 kya, until ~ 3 kya (Fig. 5). The end of cation bottlenecks, in part because vegetative this period is close to the olive domestication timeframe propagation means that perennials may not be many (~ 6.0 kya) [8] and implies a possible domestication generations removed from their ancestral genetic diver- bottleneck. Interestingly, subsequent to olive domestica- sity [54]. In order to explore this possibility, we inferred tion, an expansion of effective population size (Ne) can the demographic history of olive using SMC++ v1.15.2, a be observed. Altogether, these results suggest a mild Julca et al. BMC Biology (2020) 18:148 Page 14 of 25 bottleneck followed by a sustained population expansion (61% of the predicted Oe9 proteome) have at least one during olive domestication. SNP with nonsynonymous change, of which 10,385 are common for all the individuals (Additional file 2: Table Identification of genes under selection S8). Proteins that did not show any nonsynonymous Olive trees were domesticated for their fruits, either as a change might represent genes under particularly strong source of oil or edible fruits [56]. Genes associated with purifying selection and were enriched in several GO agronomic traits may have undergone positive selection terms (see Additional file 2: Table S9). during domestication. In order to search for genes puta- To search for signatures of selection in olive, we de- tively under positive selection in the cultivars, we first cided to compare the wild population with the two dif- classified the SNPs into intergenic, intronic, and coding. ferent groups of cultivars (cultivars_set1 and cultivars_ We further classified coding SNPs according to whether set2) using Tajima’s D and the McDonald-Kreitman they imply synonymous or nonsynonymous changes (see (MK) test for all the genes with at least 4 variable sites. “Methods”). Additional file 1: Fig. S8 shows that there is Tajima’s D had negative values for 13,258 genes for cul- a higher percentage of SNPs in intergenic regions (4.8 tivars_set1 and 7688 for cultivars_set2, of which 4464 SNPs/kb), followed by intronic (0.8 SNPs/kb) and coding are common to both groups (Additional file 2: Table regions (0.3 SNPs/kb). Moreover, the number of syn- S10). Negative Tajima’s D values indicate an excess of onymous and nonsynonymous changes is similar across low-frequency polymorphisms and can result from samples. In order to assess selection, we first measured population expansions or strong directional selection the ratio of nonsynonymous and synonymous nucleotide [59]. When we tested for functional enrichment, we ob- diversity (πN/πS) in all the individuals included in this served only two molecular functions enriched (structural study. All had similar πN/πS values, with an average of constituent of ribosome and ADP binding) in cultivars_ 0.40 (Additional file 2: Table S7), suggesting similar set1 and no enrichment for cultivars_set2. The MK test strengths of selective pressure across all sequenced indi- on all coding positions combined revealed a significant viduals. This ratio is similar to that found in other trees departure from homogeneity for both comparisons. For such as Populus nigra (0.48) [57] and Populus tricho- the first comparison, we observed that the ratio of the carpa (0.40) [58]. number of nonsynonymous to synonymous polymorphic When we analyzed the SNPs that can produce a non- sites within cultivars_set1 was significantly higher than synonymous change, including heterozygous and homo- that of the number of nonsynonymous to synonymous zygous SNPs, we found that a total of 34,060 proteins fixed sites between wild and cultivars_set1, consistent Fig. 5 SMC++ results for inferring population size histories in cultivated olives. A generation time of 20 years was used to convert coalescent scaling to calendar time Julca et al. BMC Biology (2020) 18:148 Page 15 of 25 with purifying selection (Additional file 2: Table S11). Additional file 2: Table S16), only five were found in the For the second comparison (wild vs cultivars_set2), we gene set with negative Tajima’s D, and none was present observed a significant NI < 1, consistent with positive se- in the detected selective sweep regions or the genes lection (Additional file 2: Table S11). When we per- under positive selection according to MK test. Similarly, formed the MK test in individual genes, for the a study based on transcriptomes of 68 different wild and cultivars_set1, we found 24 genes under positive selec- cultivated samples using statistical approaches (PCAdapt tion and 84 genes under negative selection. For culti- [62] and BayeScan [63]) failed to identify candidate vars_set2, we did not find genes under positive selection genes under selection associated with oil content or fruit and only one gene under negative selection (Additional size [11]. They did, however, detect ten genes as strong file 2: Table S12). candidates for selection that were associated with tran- The derived site frequency spectrum (SFS) of syn- scriptional and translational activities and to the cell onymous, nonsynonymous, fourfold degenerate, and cycle. Also, it was proposed that domestication in olive deleterious SNPs was examined for the three defined may be more related to changes in gene expression than olive populations (wild, cultivars_set1, cultivars_set2) changes in protein function in agreement with its evolu- using the subsp. cuspidata as outgroup (Additional file tionary history [11]. Further studies with more cultivated 1: Fig. S9). In all the cases, the SFS is left-shifted, with and wild samples will be needed to test this inference nonsynonymous and deleterious SNPs having higher further. density. When we compare the distributions of syn- Since some cultivars show signature of introgression onymous or fourfold degenerate SNPs with the nonsy- with sylvestris-S, we assessed whether adaptive introgres- nonymous or deleterious SNPs, we observe slightly more sion has contributed to olive domestication. When we differences in both groups of cultivars than in the wild analyzed the selective sweep regions of cultivars_set1 population, as expected from the action of purifying se- and cultivars_set2, we found that on average 42% of lection. However, these differences only reached statis- these regions overlapped introgressed regions (cultivars_ tical significance in the population with more individuals set1—43%, cultivars_set2—41%, see Additional file 2: (cultivars_set2) (p < 0.05, Mann-Whitney U test, see Table S17), with “Chemlal de Kabilye” being the cultivar Additional file 2: Table S13). with more common regions (46%). To validate our re- To further determine whether one or more proteins sults, per each cultivar, we randomly generated the same are under positive selection in cultivated olives, we number of sweep regions as the ones found in cultivars_ searched for evidence of recent selective sweeps by set1 and cultivars_set2, and check if the total length of quantifying site frequency spectrum (SFS) deviations overlap with introgressed regions is similar or higher relative to genome-wide patterns using the composite than the overlap of the predicted sweep regions. In all likelihood ratio (CLR) statistic implemented in SweeD four cultivars (“Chemlal de Kabilye,” “Megaritiki,” “Lia- [60] (see “Methods”). After merging the overlapping re- nolia Kerkyras,” and “Menya”) analyzed, we found sig- gions (see “Methods”), we kept 258 regions for cultivars_ nificant p values (p < 0.01), showing that introgressed set1 and 270 for cultivars_set2 (Additional file 2: Table regions contained selective sweeps more often than ex- S14). Shorter (< 10 kb) regions were observed in culti- pected at random. The significant overlap between re- vars_set2 (72.2% of regions are shorter than 10 kb) than gions of introgression and selective sweeps suggests a in cultivars_set1 (69.0 < 10 kb). The two groups of culti- role for adaptive introgression during olive domestica- vars have two regions larger than 50 kb, but in different tion, as it has been observed in other crops [64]. scaffolds. For cultivars_set1, 118 genes were identified in 94 sweep regions, and for cultivars_set2, 81 genes in 74 Discussion sweep regions, while the other regions were non-coding Improved assemblies for nuclear and organellar genomes in both cases. Interestingly, only 4 genes are shared by Our first reference genome assembly for Olea europaea both groups of cultivars. Among the genes found within (Oe6 version) [23] provided a needed resource for re- sweep regions, some are associated with lipid, carbohy- searchers interested in understanding the genetic basis drate, and amino acid metabolism, and others with stress of olive traits and the process of domestication. Here we tolerance, solute transport, and RNA processing (Add- improved the assembly (Oe9) by anchoring it to chro- itional file 2: Table S15). Remarkably, we did not find mosomes using a publicly available genetic map [25]. In any gene related to fatty acid metabolism and accumula- addition, we produced individual annotated assemblies tion, although this has presumably been one of the most for the plastid and mitochondrial genomes, which were important characters under selection in olive domestica- not provided as separate assemblies in earlier releases. tion. Moreover, a recent study found 19 genes associated Our comparison of our improved nuclear assembly with with five important agronomic traits in olive [61], and the recently released assembly of an oleaster from despite 14 of them being found in cv. Farga (see Turkey [24] and the identification of shared duplicated Julca et al. BMC Biology (2020) 18:148 Page 16 of 25 regions provides additional support to the previously embedded within a larger clade of other cultivars, proposed ancient polyploidization in olives [26]. Unex- whereas sister-clade relationships to other cultivars pectedly, we found apparent differences in gene content would be expected from independent domestication even using a conservative approach that is independent events. Population structure and introgression analysis of differences in annotation methodologies. These differ- discarded hypothesis (ii), as sylvestris-T individual ences could reflect shortcomings of the assemblies, showed a genetic structure similar to the nearby sampled changes occurred during domestication, or variation in cultivars and no signature of introgression. In addition, gene content among individuals, as observed in grape- the rarity of wild olive trees in the eastern Mediterra- vine, another perennial crop [65]. Further efforts in im- nean [67] together with the small number of phenotypic proving reference assemblies for the olive tree with new characters to differentiate cultivated olives from oleaster, long read technologies will definitely enhance our under- must be considered. The “Flora of Turkey” (page 156) standing on the effect of genome rearrangements, in- states that “… spontaneous seedlings may revert to var. cluding deletions of genomic regions during the sylvestris” which, despite the misuse of the term rever- evolution of olive. sion, highlights the phenotypic similarity between ferals and true oleasters in the region [67]. Earlier studies at a Reliable genetic sources: feral individuals vs. true wild large geographical scale already mentioned the difficulty trees and cultivars to find genuine oleasters in the eastern MB [51, 68]. Re- Understanding the domestication process in olives re- cent experiments confirmed this fact, analyzing putative quires careful comparison between reliably identified wild olives from Turkey [15] and Israel [12] that ap- cultivars and wild individuals. Indeed, plant choice is peared to be feral after genetic analyses. Given all these crucial before starting any genomic studies. In order to facts, we suggest that sylvestris-T is a feral individual, use reliable genetic sources, biological features of olive and we have treated it as such to avoid misleading con- germplasm have to be considered in the field: (a) high clusions. Regarding the remaining eight sylvestris in- phenotypic plasticity and resilience, (b) long reach of cluded, we hypothesize that only three are likely pure aerial transported pollen, (c) the ability to cross between sylvestris individuals, which represent only western gen- genetically distinct individuals, and (d) the long lifespan etic pool, and the other five putative oleasters are likely of the trees (often reaching several centuries and millen- introgressed individuals. Given the open debate on pos- nia). These characteristics make it extremely difficult to sible alternative centers of domestication, and the inter- tell apart pure wild individuals from those highly intro- est in tracing the genetic consequences of domestication gressed or even cultivars [66]. As stressed earlier, we in olives, the meticulous choice of additional wild popu- took great care in selecting an individual (sylvestris-S) lations, representing the eastern genetic pool, will be ne- from an isolated wild olive population, which had been cessary. Likewise, the use of not authenticated cultivars studied over the last years and meets key characteristics may lead to erroneous hypotheses. Because of this rea- that conform to a true oleaster. Our phylogenetic ana- son, we excluded the putative cultivar Dokkar [16] from lysis based on the nuclear genome places this wild indi- several analyses, given its phylogenetic position and high vidual among the earliest branching positions within the levels of introgression. This cultivar was not authenti- individuals of the subspecies europaea, and sister to cated in the original collection, and its status—wild or other wild individuals, which is fully consistent with its cultivated—is controversial, according to its morpho- ascription to a wild western Mediterranean genetic pool. logical and genetic data [9]. Thus, in order to reliably re- In stark contrast, the results for the wild sample from construct complex domestication processes, we stress Turkey [24] (sylvestris-T) were puzzling. The phylogen- the importance of careful sampling of olive trees using etic analysis placed it at a relatively shallow clade, form- ecological [69] and key morphological characters [70], ing a tight cluster with thirteen cultivars from the assisted with molecular screening of olive material, be- eastern MB. This result could be explained through four fore undertaking genome sequencing of wild and culti- hypotheses: (i) the Turkish individual identified as sylves- vated trees. tris is actually a feral olive, (ii) this tree is an oleaster highly introgressed with cultivars; (iii) considering that Genomics support for a primary domestication event the cultivated olive originated from wild oleasters, this followed by secondary events driven by introgression pattern may represent one of multiple primary domesti- from wild genetic pools cation events; or, alternatively (iv) a separate, very recent This study represents the largest phylogenetic analysis of domestication event. However, hypotheses (iii) and (iv) genome-wide sequences of the Mediterranean olives. are at odds with the observed topology in which the Our results, together with those from previous analyses well-defined clade formed by sylvestris-T and the west- [11, 16], suggest that cultivated individuals have similar ern cultivars is placed at a shallow position of the tree, nucleotide diversity as compared with wild individuals, Julca et al. BMC Biology (2020) 18:148 Page 17 of 25 being slightly higher in cultivars admixed with the west- introgression. Cultivars carrying the haplotypes E2 and ern wild genetic pool, possibly due to introgression with E3, also found in wild forms of the western MB, were local wild populations. Our demographic analyses sup- those often revealing a blueprint of introgression in our port the existence of a relatively small population size at analysis. Particularly, the cultivars harboring the E3 plas- the time of domestication, with a steady decrease in tid haplotype have a different phylogenetic history than population size preceding domestication as has been in- that of the other cultivars suggesting secondary ferred for some other crops [71, 72]. These analyses are hybridization of cultivars with wild oleasters from the consistent with earlier studies suggesting a narrow distri- Iberian Peninsula, in which a lineage similar to var. syl- bution (and hence limited population size) of oleaster vestris from Spain or Morocco acted as the maternal populations in the eastern MB over the last 150,000 line. Consistent with this, the nuclear genome of these years [8]. The demographic analysis also indicates a mild cultivars shows signs of introgression with var. sylvestris population bottleneck around 3000–14,000 years ago, from western MB. However, this introgression signal is consistent with the proposed period for primary domes- also detected in other cultivars, irrespective of the plastid tication of olives in the Levant around 6000 years ago haplotype lineage. Altogether, these results suggest that [8]. Interestingly, our analysis also suggests a rapid in- admixture with wild individuals from western popula- crease in population size following the domestication tions of var. sylvestris has been common and has taken bottleneck, likely coupled to the expansion of olive cul- place multiple times, both in the maternal and paternal tures in broader Eastern Mediterranean countries. lineages. Altogether, these results suggest that one ancestral gen- Based on plastid haplotypes of cultivars, introgression etic pool, likely deriving from a founding population from at least two different wild genetic pools, other than from eastern wild trees, is pervasive among cultivars. the one involved in the primary domestication, has con- This is consistent with common ancestry at a primary tributed to the olive germplasm. These different genetic domestication event from which all cultivars descend. pools are highly divergent from each other, as shown in However, a common primary domestication event is our phylogenetic reconstructions from plastid genomes, not incompatible with subsequent nested domestication which also suggest that the subspecies europaea might processes, perhaps driving the adaptation to local condi- have multiple origins (polyphyletic). Some other studies tions or the selection of specific traits. Considering our with plastid and nuclear markers have shown similar re- results and previous ones from numerous molecular sults [2, 3, 44]. Despite the detectable signal of introgres- techniques, the emerging pattern suggests a scenario of sion in nuclear genomes, all cultivars, including multiple events of hybridization between individuals de- introgressed ones, are monophyletic in our reconstruc- scending from this primary domestication event. These tions (considering sylvestris-T as a feral). This is also the processes have affected cultivars from the central and case of the cv. Farga, which shares a maternal ancestor western MB to different degrees. We propose that exten- with the sylvestris-S individual analyzed here. These two sive gene flow between genetically rich wild and culti- individuals are closely related when organellar genomes vated olives occurred through the expansion and are used for phylogenetic reconstruction, but cv. Farga diversification of olive crops by ancient Mediterranean appears well embedded within the cultivar clade in the cultures (Phoenicians, Romans, Arabs) [73]. It is likely nuclear genome phylogeny. This indicates that the intro- that these introgression events, followed by artificial se- gressed material in the nuclear genome has been re- lection of desired characters, resulted in the incorpor- duced through subsequent crosses following the initial ation of alleles from wild populations and facilitated the hybridization. Altogether, and in contrast to a recent creation of specific olive cultivars adapted to local study focused on olive domestication [16], the phylogen- environments. etic and admixture analyses show that rampant The origin of the introgressed genetic material can be hybridization has shaped the evolutionary history of the better inferred when the donor lineage can be traced different lineages of olives. We thus hypothesize that ol- back through the maternal line. Earlier studies have ives represent a domestication continuum [6]. Following shown the presence of clearly distinct haplotype groups a primary domestication event in the eastern MB, there (E1, E2, and E3) among cultivars [17, 44]. E1 is present were additional hybridizations with local wild genetic in wild trees from the eastern MB and it is likely the sig- pools throughout the MB. Although there may be some nature of the primary domestication event, since it is evidence and rationale for calling these hybridizations shared by the 90% of the current olive cultivars [18]. secondary domestication events [12], all of the germ- This is particularly true for individuals sampled from re- plasm in our analyses contain some remnant of the pri- gions close to the origin of domestication, cultivars from mary event. We note, however, increased genomic Syria, Iran, Turkey, and Lebanon and the feral individual sampling, particularly of truly wild populations from the from Turkey, which show little or no signs of Levant, is needed to help describe these complex Julca et al. BMC Biology (2020) 18:148 Page 18 of 25 patterns of evolution in Olea europaea across the nu- Conclusions merous areas of the Mediterranean basin. The comparison of the improved reference genome of Olea europaea cv. “Farga” (version Oe9) with the other two available genomes and all sequenced individuals Lack of a clear domestication syndrome showed a different gene content among individuals. Cultivated olives have undergone a complex domestica- Phylogenetic and introgression analysis helped to un- tion process, which has led to morphological and physio- cover the feral origin of six samples. Remarkably, our re- logical changes. The main observable traits selected sults suggest that the available genome of var. sylvestris during the transition between oleasters and cultivars are corresponds to a feral individual. fruit size and oil content [13, 56, 70]. The domestication Genes positively selected in cultivated olives were as- scenario described in the previous section, which is sociated mainly with response to biotic and abiotic punctuated by hybridization, may make it difficult to de- stresses. We did not find evidence of selection of genes tect genes selected during the process of domestication. associated with fruit size and oil content. We used three different approaches (McDonald-Kreit- Rampant hybridization has shaped the evolutionary man test, Tajima’s D, and selective sweeps prediction) to history of the different lineages of olives. After the pri- detect sequence changes likely associated with selection. mary domestication in the eastern Mediterranean basin, We detected different sets of genes putatively selected there were numerous secondary events across most for cultivars_set1 and cultivars_set2 with few genes over- countries of southern Europe and northern Africa, often lapping the different methods. In general, genes posi- involving genetic admixture with genetically rich wild tively selected in cultivated olives were associated with a populations, particularly from the western Mediterra- response to biotic and abiotic stresses. Interestingly, five nean Basin. genes detected with negative Tajima’s D are related to genes recently described as important for fruit weight Methods and stone weight [61]. However, further analyses are Scaffolding of the cultivated olive genome using a needed to ascertain whether they indeed play a role in linkage map trait selection during domestication (see [74] for second- A new, improved version of the O. europaea genome as- ary compounds). Previous analyses detected few signs of sembly (Oe9) was produced by anchoring the Oe6 ver- positive selection affecting protein-coding regions, but sion [23] to a publicly available genetic map [25] using they did detect differences in expression levels between ALLMAPS [75]. First, we took the intersection of 7042 cultivars and wild individuals for specific genes [11]. markers for which a sequence was provided and mapped This led to the conclusion that selection may have acted them with BWA v0.7.15 [76] to the Oe6 olive genome on non-coding regions that drive gene expression. How- assembly. Filtering for minimum mapping quality 20 and ever, given the difficulty of controlling for other factors fewer than 10 mismatches, we obtained 5780 mappings. driving gene expression (different periods of the year, Intersecting these mapped markers with the 3404 local environmental conditions), we believe this result markers placed in the genetic map resulted in a set of must be viewed with caution. Also, a recent study 2759 markers with both a genetic map position and an showed that the activation of transposable elements near unambiguous physical location in the Oe6 assembly. genes of important agronomic traits may play a role in ALLMAPS was then run with default parameters. A total olive domestication [16]. Alternative approaches are of 2362 markers were considered unique, of which 2134 needed to detect alleles whose presence in the cultivars were anchored and 228 unplaced. were selected through domestication. We cannot rely on models that assume only vertical evolution, but Functional annotation of the cultivated olive genome would rather search for shared conserved or intro- In order to give more functional insight to our analysis, gressed regions across different cultivars sharing simi- we decided to improve the functional annotation of the lar phenotypes [70]. To be effective, a much larger genes present in the Oe6 version of the genome by run- sampling of genomic sequence will be needed for ning Blast2go [77], which in turn ran a BLASTP [78] such an approach. In addition, assemblies from repre- search against the non-redundant database (April 2018) sentative lineages of reliable wild olives will help to and Interproscan [79] to detect protein domains on the better trace the origin of different introgressed re- annotated proteins. gions in the genomes of cultivars. Overall, our new phylogenomic and genetic analyses of whole genome Comparison of the europaea and sylvestris genomes sequences show evidence for a complex process of re- To compare the gene completeness of the four available ticulation disrupting historical isolation in the course genome assemblies of O. europaea subsp. europaea (two of olive domestication. var. europaea of cv. Farga, one cv. Picual, and one var. Julca et al. BMC Biology (2020) 18:148 Page 19 of 25 sylvestris), the BUSCOv3 (Benchmarking Universal they help to provide a temporal polarity to the data (e.g., Single-Copy Orthologs) program was run using the distinguish between ancestral and derived characters). embryophita_odb9 database made of 1440 conserved These samples broadly covered the geographical distri- orthologous genes in plants (Additional file 2: Table S1). bution of the species in the MB (see Table 1). The au- To better characterize gene completeness differences be- thenticity of these cultivars was previously substantiated tween the cvs. Farga and Picual genomes assemblies that through molecular and morphological markers [9]. The had similar BUSCO results, two RNAseq datasets ob- sequenced oleaster (referred to as sylvestris-S from now tained from immature cv. Farga olives (ERR1406351) onwards) was collected in the North of Spain. All the and cv. Arbequina seeds (SRR6003535) were aligned to samples were collected according to the local, national, the two assemblies with STAR. Total percentage of reads or international guidelines and legislation. mapped, as well as percentage of unique mappings, has The DNA of all individuals was extracted as described been reported (Additional file 2: Table S2). in [23], and their genomes were sequenced using Illu- Additionally, to compare the three available genome mina HiSeq 2000 paired-end technology to a sequencing assemblies of var. europaea (“Farga” v. Oe6 and Oe9, depth ranging from 24 to 34× at the CNAG-CRG se- “Picual”) and sylvestris, we plotted their cumulative as- quencing facilities, as described for the reference gen- sembly lengths using the python package seaborn [80]. ome [23]. In addition to these individuals, we used The predicted protein-coding gene sequences of the two publicly available data of the reference genome sequence assemblies were compared using a search with BLASTN of the olive cv. Farga [23]. “Farga” is a cultivar from [27]. Results were filtered using cutoffs of 80% identity Catalonia (eastern Spain) with the E3.1 plastid haplotype and e-value < 1e−5. For the genes that did not have a and previously classified as representative of the central hit, we analyzed whether they were covered by reads of MB cultivated genepool (Table 1). We also used a pub- the other sequenced accession. For this step, we first lished assembly of O. europaea var. sylvestris (referred as mapped the reads of each genome against the other sylvestris-T from now onwards) [24], and 49 samples using BWA 0.7.6a-r433 [76]. We considered a gene as from a recent publication [16], and downloaded fourteen individual-specific (e.g., europaea-specific) if it did not plastid genomes of the Mediterranean olive from the pass a filter of coverage > 20 over 50% of the gene length NCBI database (see Table 1). in the other genome (e.g., sylvestris). With the aim of achieving a more conservative set of individual-specific genes, we applied a more stringent filter of coverage (> Organelle assemblies 5). For the list of detected individual-specific genes The available reference genome sequence does not in- found in europaea, we mapped the reads of all the var- clude separate scaffolds for mitochondrial and plastid ieties and an outgroup (subsp. cuspidata) using the same genomes [23]. Here, we assembled and annotated both approach. Bedtools v2.26.0 was then used to detect organellar genomes of the cv. Farga using paired-end coverage and a gene was assumed to be present in a var- (PE) and mate-pair (MP) data from the reference gen- iety when at least 90% of the coding region was mapped ome sequence project [23]. Briefly, all genome shotgun by reads (coverage filter > 5). The set of genes with no Illumina reads were mapped using BWA v0.7.13-r1126 reads mapping from any other variety was scanned for [76] to a reference plastid sequence (NC_013707), and a additional contamination. In total, 104 genes were dis- mitochondrial sequence (MG372119). Then reads carded due to the fact that they did only have hits in were filtered allowing only those that mapped in other non-plant species. proper pairs with a hard and soft clipping for a max- Finally, in order to search for the functions of the imum of the 25% of the total length of the read. The genes without a hit in the other genome (individual-spe- plastid genome was assembled using NOVOPlasty cific), we performed a BLAST search against the NCBI v2.6.3 [81]. The mitochondrial genome was first as- non-redundant database, and the same cutoffs as de- sembled with SPAdes v3.10.0 [82]. Then, this initial scribed before. assembly was scaffolded using SSPACE_Basic v2.0 [83] using PE and MP libraries. Finally, gaps were Genome sequences filled with GapFiller v1–10 [84]. We sampled and sequenced twelve genomes of O. euro- The plastid and mitochondrial genomes were annotated paea: ten cultivars (“Arbequina”, “Beladi”, “Picual,” “Sor- by BLAST searches against previously annotated plastid ani,” “Chemlal de Kabilye,” “Megaritiki,” “Lechin de (FN996972, MG255763) and mitochondrial (MG372119, Sevilla,” “Lechin de Granada,” “Frantoio,” and “Koro- KX545367) genomes. Gene structures (i.e., intron–exon neiki”), one oleaster (sylvestris-S), and one subsp. cuspi- boundaries) were defined using Exonerate v2.47.3 [85] data to be used as an outgroup. The inclusion of using the “protein2genome” model. Annotations of tRNA outgroups in phylogenomic analyses is recommended as genes were verified using tRNAscan-SE [86]. Julca et al. BMC Biology (2020) 18:148 Page 20 of 25 Detection of single nucleotide variants in a pseudo-alignment of all the sequenced genomes. To assess nucleotide diversity across sequences of the Specifically, for the nuclear dataset, we included only Mediterranean olive at the nuclear, plastid, and mito- homozygous SNPs. For the plastid genomes, we included chondrial levels, we called single-nucleotide polymor- additional sequences by aligning our genomes with the phisms (SNPs) using the cv. Farga genome as a reference database genomes (see Table 1) using MAFFT v7.305b for all the cases (for the nuclear genome we used the [93]. All these alignments were trimmed using trimAl Oe9 version). Sequenced reads from each individual v1.4 [94] with the options -st 1 and -complementary, in were mapped against the reference genome using BWA order to remove all the non-informative positions. The 0.7.6a-r433 [76], and SNPs were identified with GATK final alignment had 12,862,844 variable positions for the HaplotypeCaller v4.0.8.1 [87], setting ploidy to 2, and nuclear genome, 327 for the plastid genome, and 4168 using thresholds for mapping quality (MQ > 40), quality for the mitochondrial genome. Phylogenetic trees were by depth (QD > 2), Fisher strand bias (FS < 60), mapping reconstructed from these alignments using RAxML quality rank sum test (MQRankSum > − 12.5), read pos v8.1.17 [95] and the GTR model as it is the most fre- rank sum test (ReadPosRankSum > − 8), strand odds ra- quent evolutionary model found in previous studies. tio (SOR < 3), read depth of coverage (DP ≥ 10) in at Support values were calculated based on 100 bootstrap least the 80% of individuals, and allelic depth (AD ≥ 5) in searches using the rapid bootstrapping as implemented at least the 80% of individuals. Sites with missing alleles in RAxML. Additionally, for the nuclear data, we recon- and spanning a deletion were also removed. Finally, structed a phylogenetic network using SplitsTree4 VCFtools v0.1.17 [88] was used to filter out positions ac- v4.14.5 and the NeighborNet approach [96]. Since we cording to the number of alleles (--max-alleles 2) and had six cultivars (“Koroneiki,” “Lechin de Sevilla,” “Arbe- minor allele frequency (--maf 0.008). quina,” “Chemlal de Kabylie,” “Picual,” “Frantoio”) that had two samples and in the nuclear tree cluster together, we decided to only keep the individuals that were se- Admixture mapping quenced in this project. After this filter, we keep 46 cul- We used Bayesian hierarchical clustering and principal tivars, 10 sylvestris, and 1 subsp. cuspidata, which were component analysis (PCA) of genetic variance to identify used for all the analysis. population structure without a priori grouping assump- tions among individuals of the subsp. europaea. For the Bayesian hierarchical clustering, we used the Structure Analysis of introgression with SNP data software v2.3.4 (89) and because of the large number of The ABBA-BABA test [49] was used to search for evi- polymorphic positions in the nuclear genomes of the O. dence of introgression among the olive samples. Dsuite europaea samples, and computational limitations, we v0.1 r3 [50] was employed to calculate the D-statistic generated ten subsets of 100,000 randomly chosen poly- from nuclear SNP data for all subsets of quartets that morphic positions without overlaps and analyzed them were compatible with the previously reconstructed in parallel. Structure was run with 100,000 generations phylogenetic tree (see above). Each quartet includes the of “burn-in” and 100,000 Markov chain Monte Carlo tree subsp. europaea individuals plus the subsp. cuspi- (MCMC) iterations after burn-in for increasing K values data as outgroup (Additional file 2: Table S18). For mul- ranging from 1 to 6, considering independent alleles and tiple hypotheses testing, we applied a false discovery rate admixture of individuals. Simulations were repeated five correction to the p values [97]. Then a heatmap showing times for each value of K. The optimal number of gen- the D-statistic and its p value was plotted for all pairs of etic clusters was determined using the ΔK method [89] individuals using the plot_d.rb script (https://github. using the software Structure Harvester [90]. The optimal com/mmatschiner/tutorials/blob/master/analysis_of_ K value was visualized with DISTRUCT v1.1 [91]. For introgression_with_snp_data/src/plot_d.rb). The Dsuite PCA analyses, we used PLINK v2.00a2.3LM [92] and the Dinvestigate tool was used in sliding windows of 5000 complete set of SNPs (24,724,756 positions). Finally, SNPs, incremented by 1000 SNPs to test whether the D- population differentiation (weighted FST) was calculated statistic is homogeneous or variable throughout the gen- using VCFtools v0.1.17 [88]. ome. It was tested in the trios with the strongest signals of introgression: “Beladi”-“Dokkar”-sylvestris-S, “Beladi”- Phylogenetic analysis “Chemlal de Kabilye”-sylvestris-S, “Beladi”-”Megaritiki”- Phylogenetic trees were reconstructed using SNP data sylvestris-S, “Beladi”-”Menya”-sylvestris-S, “Beladi”-”Lia- from nuclear, plastid, and mitochondrial genomes separ- nolia Kerkyras”-sylvestris-S. Then, the overlapping re- ately. In each case, the genome sequence of each se- gions with D-statistic > 0.5 were joined per each of the quenced individual was obtained by replacing the SNP five cultivars analyzed and genes present in these regions positions in the respective reference genome, resulting were used for enrichment analysis. Julca et al. BMC Biology (2020) 18:148 Page 21 of 25 Nucleotide diversity population (see section “Nucleotide diversity”). To evalu- Pairwise nucleotide diversity (π) was calculated using ate if the ratio of synonymous to nonsynonymous poly- VCFtools v0.1.17 [88] per window of 20 kb. The samples morphisms within individuals of cultivars (cultivars_set1 were grouped together into three datasets: oleaster (only and cultivars_set2 independently) will be similar to the ra- the sylvestris with 100% of cluster 1 as predicted with tio of synonymous to nonsynonymous divergences be- Structure (Additional file 2: Table S5)), cultivars_set1 (cul- tween cultivars and wild populations (fixed differences), tivars with genetic cluster 2), cultivars_set2 (cultivars with we performed a McDonald-Kreitman test for all coding clusters 1 and 2). The sylvestris with mixed clusters regions [102]. We also reported the neutrality index (NI), (W2R74, W9R302, W1R198, W4R183, W8R225) and syl- which shows the directionality of the McDonald-Kreitman vestris-T were excluded because of their possible feral ori- test, and the p value was calculated using Fisher’s exact gin. The cv. Dokkar was also excluded for its phylogenetic test, which evaluates if the differences of the ratios are sig- position and high levels of introgression (see Fig. 2a,b). nificant. An NI value > 1 is consistent with negative selec- tion, while an NI value < 1 is consistent with positive Demographic history of cultivated olives selection. For multiple hypotheses testing, we applied a To estimate population size histories in olive, we false discovery rate correction to the p values [97]. How- employed SMC++ v1.15.2 [55], which is capable of analyz- ever, after the correction, none of the genes shows signifi- ing unphased genomes. The dataset included only 45 cul- cant p values, but anyway we decided to show the few tivars, we decided to exclude “Dokkar” for its phylogenetic genes with p < 0.05. We also generated the derived site fre- positions and high levels of introgression (see Fig. 2a,b). quency spectrum (SFS) for synonymous, nonsynonymous, First, we masked all regions larger than 5 kb with a cover- fourfold degenerate, and deleterious sites using subsp. cus- age < 10× in at least one of the samples included using pidata to determine ancestral state. Then we compared bedtools v2.26.0 [98]. Then, we ran SMC++ using default the SFS distribution of probable neutral evolving sites parameters and setting the T1, the most recent time point (synonymous and fourfold degenerate) with no neutral for population size history inference, to 10. Finally, a gen- evolving sites (nonsynonymous sites and deleterious) eration time of 20 years [12] and a mutation rate of 7.77 × using the Mann-Whitney U test. Additionally, an ap- 10− 9 mutations per nucleotide per generation [99] were proach based on the site frequency spectrum (SFS). used to convert the scaled times and population sizes into SweeD v3.1 [60] was used to identify selective sweeps in real times and sizes. olive. This program is based on Sweepfinder [103] and uses a composite likelihood ratio (CLR) to identify loci SNP characterization showing a strong deviation in the site frequency spectrum Nuclear SNPs were classified according to their genomic (SFS) toward rare variants. We used the subsp. cuspidata location as intergenic, intronic, and coding. Coding as an outgroup to infer ancestral alleles. SweeD was run SNPs were further classified into synonymous and non- separately for each scaffold and grid as the only parameter. synonymous, according to the implied change in the re- The grid parameter was calculated per scaffold in order to spective codon. We also differentiate the fourfold have a measure of the CLR every 5 kb (size of the scaf- degenerate positions and deleterious SNPs (if a nonsy- fold/5000). Outliers were defined as the 0.5% with the nonymous SNP produces a change in the start or stop most extreme p values. Closer regions, with less than 10 codon). For the heterozygous positions, if at least one of bp distance, were collapsed. Finally, for all the protein lists the alleles was nonsynonymous, we classified the pos- generated, we performed a GO term enrichment analysis ition as nonsynonymous. GO term enrichment analyses using FatiGO [100]. of the proteins without nonsynonymous SNPs was cal- culated using FatiGO [100]. To investigate the variation Identification of overlapping sweep and introgressed of nonsynonymous and synonymous SNPs within coding regions regions, we compared nonsynonymous changes per non- In order to test whether adaptive introgression played an synonymous site (πN) to synonymous changes per syn- important role in olive domestication, we used the sweep onymous site (πS) using a synonymous: nonsynonymous regions generated for cultivars_set1 and cultivars_set2 site ratio of 1:3 [101]. (Additional file 2: Table S14) and compared with the predicted introgressed regions for cultivars “Chemlal de Identification of genes under selection Kabilye”, “Megaritiki”, “Lianolia Kerkyras” and “Menya” To detect protein-coding genes that have potentially (Additional file 2: Table S6). From the sweep regions, undergone selection among the cultivated individuals, we only those with 100% overlap were selected. To validate used different approaches. Tajima’s D neutrality test was our results, we calculated an empirical p value for each performed on coding sequences of all genes with at least 4 cultivar that tests whether random generated regions SNPs using VCFtools v0.1.17 [88] for each defined olive have larger overlap with the introgressed regions than Julca et al. BMC Biology (2020) 18:148 Page 22 of 25 the predicted sweep regions. Briefly, the same number of of the var. sylvestris from Turkey. sweep regions were randomly generated with similar Table S2. Total percentage of reads mapped from two RNAseq libraries, as sizes, and then we check if the total length of overlap well as, percentage of unique mappings to the Oe9 and picual genomes. Table S3. List of unique genes of europaea and sylvestris with their with introgressed regions (sampled) is the same or larger homologous function. Table S4. General characteristics of the plastid and than that of the sweep regions (observed). We repeated mitochondrial genomes of the cultivar Farga. Table S5. Admixture this procedure 1000 times and calculated the p value by coefficient (Q) of each individual per cluster. This table was used to create the Fig. 2b. Table S6. Introgressed regions per each cultivar. Columns in dividing the number of times where the length of sam- order show: cultivar name (‘Chemlal de Kabilye’, ‘Megaritiki’, ‘Lianolia pled regions ≥ observed regions by the number of Kerkyras’, ‘Menya’, and ‘Dokkar’), scaffold, start of the region, end of the repetitions. region, D-statistic, genes present in the region. Table S7. Number of syn- onymous, nonsynonymous, four fold degenerate, affecting stop codons, af- fecting start codons SNPs for homozygous and heterozygous positions per individual. TheπN/πS ratio of homozygous SNPs, πN/πS ratio of heterozy- Supplementary information gous SNPs, and πN/πS ratio of total number of SNPs is also shown. Table Supplementary information accompanies this paper at https://doi.org/10. S8. Number of proteins with SNPs and with nonsynonymous SNPs per indi- 1186/s12915-020-00881-6. vidual. Table S9. GO terms enriched in the list of proteins that do not have a nonsynonymous SNP. First column shows the term category, the second, Additional file 1 Fig. S1. Genome comparison of sylvestris and the GO term, the third, the term level, the fourth, the p-value, and the fifth, europaea. (a) Cumulative genome length per scaffold ranked in order of the term name. Table S10. List of genes with negative Tajima’s D. The col- size for the genome assembly of europaea (Oe9 - green, Oe6 - olive, umns in order indicate: the group of cultivars (cultivars_set1 or cultivars_- picual - orange) and sylvestris (blue). A straight vertical line represents a set2), gene, scaffold, number of SNPs, Tajima’s D. Table S11. Result of perfect genome assembly. The horizontal plateaus indicate many small McDonald and Kreitman test for all coding regions of wild vs cultivars_set1 scaffolds. The top right end of each curve shows the total number of and wild vs cultivars_set2. Table S12. Result of McDonald and Kreitman scaffolds. (b) Syntenic plot of the genome of europaea against sylvestris test for each gene. The columns in order show: the comparison for the test generated by SynMap. (c) Histogram of log 10 transformed Ks values of (wild vs cultivars_set1 or wild vs cultivars_set2), number of fixed nonsynon- syntenic gene pairs identified as calculated by SynMap. Fig. S2. Plastid ymous sites (Fn), number of fixed synonymous sites (Fs), number of poly- genome of the cultivar Farga. Protein coding genes are shown in green, morphic nonsynonymous sites (Pn), number of polymorphic synonymous rRNAs in light blue, and tRNAs in purple. The SNPs are shown per each sites (Ps), the ratio of fixed nonsynonymous sites per synonymous sites (Fn/ individual sequenced in this study in the following order starting from Fs), the ratio of polymorphic nonsynonymous sites per synonymous sites outside: ‘Arbequina’, ‘Picual’, ‘Beladi’, ‘Sorani’, ‘Koroneiki’, ‘Frantoio’, ‘Lechin (Pn/Ps), the neutrality index (NI), p-value (Fisher’s exact tests), homologous de Granada’, ‘Lechin de Sevilla’, ‘Megaritiki’, ‘Chemlal de Kabilye’, sylvestris- function. Table S13. Derived site frequency of synonymous, nonsynon- T, sylvestris-S. Fig. S3. Mitochondrial genome of the cultivar Farga. ymous, deleterious, and four fold degenerate sites. The p-values using the Protein coding genes are shown in green, rRNAs in light blue, and tRNAs Mann-Whitney U test of the comparison of neutral and no neutral sites are in purple. The SNPs are shown per each individual sequenced in this also shown. Table S14. Selective sweeps for cultivars_set1 and cultivars_- study in the following order starting from outside: ‘Arbequina’, ‘Picual’, set2. The genes that are present in the regions are indicated. Table S15. ‘Beladi’, ‘Sorani’, ‘Koroneiki’, ‘Frantoio’, ‘Lechin de Granada’, ‘Lechin de List of proteins in regions with selected sweeps and their associated func- Sevilla’, ‘Megaritiki’, ‘Chemlal de Kabilye’, sylvestris-T, sylvestris-S. Fig. S4. tion. Table S16. Blast results of the 19 genes of sylvestris-T against cv. Farga. Homozygous SNP distribution along the nuclear genome. The SNPs are The results were filtered by %identity > 90 and e-value<1e-5. Asterisk mark shown in windows of 100 Kb. Cultivars are plotted in green and sylvestris genes with negative values of Tajima’s D in cultivars_set1 (*), in cultivars_- in blue. Since cv. Farga was used as a reference genome, we do not set2 (**), and in both sets of cultivars (***). Table S17. Selective sweeps expect homozygous SNPs for this sample. Fig. S5. Heterozygous SNP present in introgressed regions of ‘Chemlal de Kabilye’, ‘Megaritiki’, ‘Lianolia distribution along the nuclear genome. The SNPs are shown in windows Kerkyras’, and ‘Menya’. Table S18. D-statistic of all trios of subsp. europaea of 100 Kb. Cultivars are plotted in green and sylvestris in blue. Fig. S6. analysed. The columns show the name of the individual 1, individual 2, Principal component analysis (PCA) based on 24,724,756 nuclear SNPs. individual 3, D-statistic, p-value, and adjusted p-value. The first two principal components (PC1 and PC2) are plotted, which explained 35.17% and 11% of the total variance, respectively. The individuals are shown with three different colors according to the results Acknowledgements of structure (Fig. 2b): blue for individuals composed mainly of cluster 1, We want to acknowledge Guillaume Besnard, Manuel Nogales, Manuel green for individuals composed mainly of cluster 2, and orange for Sánchez, and Carlos Garcia-Verdugo for providing some of the plant samples mixed individuals. The cultivars are identified by a circle and the oleaster used in this study. The authors want to thank Danelle Seymour and Rebecca by a cross (x). Fig. S7. Patterns of introgression. (a) D statistics across the Gaout for preparing some of the libraries. TG and PV acknowledge support largest scaffold (Oe9_LG02) in windows of 5Kb and step size of 1Kb for from Banco Santander for the Olive Genome Sequencing project. the analysis of the trios: ‘Beladi’-‘Dokkar’-sylvestris-S, ‘Beladi’-‘Chemlal de Kabilye’-sylvestris-S’, ‘Beladi’-‘Megaritiki’-sylvestris-S, ‘Beladi’-‘Menya’-sylvestris- Authors’ contributions S, ‘Beladi’-‘Lianolia Kerkyras’-sylvestris-S. The horizontal red line marks D = IJ, MMH, FC, JGG, and TSA performed bioinformatics analysis. IJ and TG 0.5. (b) Genomic regions introgressed with a D > 0.5 in ‘Megaritiki’ (grey), drafted a first version of the manuscript to which all authors contributed. IJ, ‘Menya’ (blue), ‘Lianolia Kerkyras’ (orange), ‘Dokkar’ (green), and ‘Chemlal MMH, FC, JGG, BSG, CMD, IGG, TSA, PV, and TG read and approved the final de Kabilye’ (red). Percentages below the cultivar names represent the version of the manuscript. TSA and IGG supervised the sequencing, percentage of the genomic region that is introgressed in the cultivar. assembly, and annotation. CMD, BSG, and PV provided samples, designed Fig. S8. Number of homozygous and heterozygous SNPs (SNPs/Kb) in the selection of individuals, and provided data and important scientific the intergenic, intronic and coding region of the genome. The coding insights in the interpretation of the results. TG supervised the study. The region was divided according to the changes that the allele can produce author(s) read and approved the final manuscript. (synonymous and nonsynonymous). Fig. S9. Derived site frequency spectrum (SFS) of synonymous, four fold degenerate, non-synonymous, and deleterious SNPs. (a) Wild populations, (b) Cultivars_set1, (c) Culti- Funding vars_set2. Subspecies cuspidata was used as an outgroup for identifying This research has received funding from the European Research Council the derived allele. (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 724173; RETVOLUTION)”. IJ was supported in Additional file 2 Table S1. General characteristics of the four genome assemblies of the cv. Farga (Oe6, Oe9), picual, and the genome assembly part by a grant from the Peruvian Ministry of Education: “Beca Presidente de la República” (2013-III). Julca et al. BMC Biology (2020) 18:148 Page 23 of 25 Availability of data and materials using SSR and morphological markers. Tree Genet Genomes. 2014;10:141– All data generated or analyzed during this study are included in this 55. https://doi.org/10.1007/s11295-013-0671-3. published article, its supplementary information files and publicly available 10. Besnard G, Rubio de Casas R. Single vs multiple independent olive repositories. The new version of the O. europaea var. europaea genome and domestications: the jury is (still) out. New Phytol. 2016;209:466–70. https:// the resequencing data produced in this project have been deposited in the doi.org/10.1111/nph.13518. ENA (European Nucleotide Archive) under the accession number 11. Gros-Balthazard M, Besnard G, Sarah G, Holtz Y, Leclercq J, Santoni S, et al. PRJEB35540. The genome assembly and annotation are also available in Evolutionary transcriptomics reveals the origins of olives and the genomic https://denovo.cnag.cat/olive, together with a JBrowse genome browser and changes associated with their domestication. Plant J. 2019; https://doi.org/ a Blast Server. 10.1111/tpj.14435. 12. Diez CM, Trujillo I, Martinez-Urdiroz N, Barranco D, Rallo L, Marfil P, et al. Ethics approval and consent to participate Olive domestication and diversification in the Mediterranean Basin. New Not applicable. Phytol. 2015;206:436–47. https://doi.org/10.1111/nph.13181. 13. Breton C, Terral J-F, Pinatel C, Médail F, Bonhomme F, Bervillé A. The origins of the domestication of the olive tree. C R Biol. 2009;332:1059–64. https:// Consent for publication doi.org/10.1016/j.crvi.2009.08.001. Not applicable. 14. Díez CM, Gaut BS. The jury may be out, but it is important that it deliberates: a response to Besnard and Rubio de Casas about olive domestication. New Competing interests Phytol. 2016;209:471–3. https://doi.org/10.1111/nph.13780. The authors declare that they have no competing interests. 15. Yoruk B, Taskin V. Genetic diversity and relationships of wild and cultivated olives in Turkey. Plant Syst Evol. 2014;300:1247–58. https://doi.org/10.1007/ Author details s00606-014-1002-3. 1 Centre for Genomic Regulation (CRG), The Barcelona Institute of Science 16. Jiménez-Ruiz J, Ramírez-Tejero JA, Fernández-Pozo N, Leyva-Pérez M de la and Technology, Dr. Aiguader 88, 08003 Barcelona, Spain. 2Department of O, Yan H, Rosa R de la, et al. Transposon activation is a major driver in the Experimental and Health Sciences, Universitat Pompeu Fabra (UPF), 08003 genome evolution of cultivated olive trees ( Olea europaea L.). Plant Barcelona, Spain. 3Universitat Autònoma de Barcelona (UAB), 08193 Genome. 2020;:e20010. doi:https://doi.org/10.1002/tpg2.20010. Barcelona, Spain. 4Present address: Barcelona Supercomputing Centre 17. Besnard G, Hernández P, Khadari B, Dorado G, Savolainen V. Genomic (BSC-CNS), and Institute for Research in Biomedicine (IRB), Barcelona, Spain. profiling of plastid DNA variation in the Mediterranean olive tree. BMC Plant 5 CNAG-CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Biol. 2011;11:80. https://doi.org/10.1186/1471-2229-11-80. Science and Technology (BIST), Baldiri i Reixac 4, 08028 Barcelona, Spain. 18. Besnard G, El Bakkali A, Haouane H, Baali-Cherif D, Moukhli A, Khadari B. 6 Department Ecology and Evolutionary Biology, University of California Irvine, Population genetics of Mediterranean and Saharan olives: geographic Irvine, CA 92697, USA. 7Agronomy Department, University of Cordoba, 14071 patterns of differentiation and evidence for early generations of admixture. Cordoba, Spain. 8Royal Botanical Garden of Madrid. Consejo Superior de Ann Bot. 2013;112:1293–302. https://doi.org/10.1093/aob/mct196. Investigaciones Científicas (CSIC), 28014 Madrid, Spain. 9Institució Catalana de 19. Haouane H, El Bakkali A, Moukhli A, Tollon C, Santoni S, Oukabli A, et al. Recerca i Estudis Avançats (ICREA), Pg. Lluís Companys 23, 08010 Barcelona, Genetic structure and core collection of the World Olive Germplasm Bank Spain. of Marrakech: towards the optimised management and use of Mediterranean olive genetic resources. Genetica. 2011;139:1083–94. https:// Received: 28 March 2020 Accepted: 27 September 2020 doi.org/10.1007/s10709-011-9608-7. 20. Díez CM, Imperato A, Rallo L, Barranco D, Trujillo I. Worldwide core collection of olive cultivars based on simple sequence repeat and morphological References markers. Crop Sci. 2012;52:211. https://doi.org/10.2135/cropsci2011.02.0110. 1. Green PS. A revision of olea L. (oleaceae). Kew Bull. 2002;57:91. https://doi. 21. El Bakkali A, Haouane H, Moukhli A, Costes E, Van Damme P, Khadari B. org/10.2307/4110824. Construction of core collections suitable for association mapping to 2. Besnard G, Rubio de Casas R, Christin P-A, Vargas P. Phylogenetics of Olea optimize use of Mediterranean olive (Olea europaea L.) genetic resources. (Oleaceae) based on plastid and nuclear ribosomal DNA sequences: tertiary PLoS One. 2013;8:e61265. https://doi.org/10.1371/journal.pone.0061265. climatic shifts and lineage differentiation times. Ann Bot. 2009;104:143–60. 22. Besnard G, Baradat P, Breton C, Khadari B, Bervillé A. Olive domestication https://doi.org/10.1093/aob/mcp105. from structure of oleasters and cultivars using nuclear RAPDs and 3. Vargas P, Kadereit JW. Molecular fingerprinting evidence (ISSR, Inter-Simple mitochondrial RFLPs. Genet Sel Evol. 2001;33:S251–68. https://doi.org/10. Sequence Repeats) for a wild status of Olea europaea L. (Oleaceae) in the 1186/BF03500883. Eurosiberian North of the Iberian Peninsula. Flora. 2001;196:142–52. https:// 23. Cruz F, Julca I, Gómez-Garrido J, Loska D, Marcet-Houben M, Cano E, et al. doi.org/10.1016/S0367-2530(17)30029-4. Genome sequence of the olive tree, Olea europaea. Gigascience. 2016;5:29. 4. Rubio de Casas R, Besnard G, Schönswetter P, Balaguer L, Vargas P. https://doi.org/10.1186/s13742-016-0134-5. Extensive gene flow blurs phylogeographic but not phylogenetic signal in 24. Unver T, Wu Z, Sterck L, Turktas M, Lohaus R, Li Z, et al. Genome of wild Olea europaea L. Theor Appl Genet. 2006;113:575–83. https://doi.org/10. olive and the evolution of oil biosynthesis. Proc Natl Acad Sci U S A. 2017; 1007/s00122-006-0306-2. 114:E9413–22. https://doi.org/10.1073/pnas.1708621114. 5. Zohary D, Spiegel-Roy P. Beginnings of fruit growing in the old world. 25. İpek A, İpek M, Ercişli S, Tangu NA. Transcriptome-based SNP discovery by Science. 1975;187:319–27. https://doi.org/10.1126/science.187.4174.319. GBS and the construction of a genetic map for olive. Funct Integr 6. Kaniewski D, Van Campo E, Boiy T, Terral J-F, Khadari B, Besnard G. Primary Genomics. 2017;17:493–501. https://doi.org/10.1007/s10142-017-0552-1. domestication and early uses of the emblematic olive tree: palaeobotanical, 26. Julca I, Marcet-Houben M, Vargas P, Gabaldón T. Phylogenomics of the olive historical and molecular evidence from the Middle East. Biol Rev Camb tree (Olea europaea) reveals the relative contribution of ancient allo- and Philos Soc. 2012;87:885–99. https://doi.org/10.1111/j.1469-185X.2012.00229.x. autopolyploidization events. BMC Biol. 2018;16:15. https://doi.org/10.1186/ 7. Terral J-F, Alonso N, RB i C, Chatti N, Fabre L, Fiorentino G, et al. Historical s12915-018-0482-y. biogeography of olive domestication (Olea europaea L.) as revealed by 27. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. geometrical morphometry applied to biological and archaeological material. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10:421. J Biogeogr. 2004;31:63–77. https://doi.org/10.1046/j.0305-0270.2003.01019.x. https://doi.org/10.1186/1471-2105-10-421. 8. Besnard G, Khadari B, Navascués M, Fernández-Mazuecos M, El Bakkali A, 28. de Abreu-Neto JB, Turchetto-Zolet AC, de Oliveira LFV, Zanettini MHB, Arrigo N, et al. The complex history of the olive tree: from Late Quaternary Margis-Pinheiro M. Heavy metal-associated isoprenylated plant protein diversification of Mediterranean lineages to primary domestication in the (HIPP): characterization of a family of proteins exclusive to plants. FEBS J. northern Levant. Proc Biol Sci. 2013;280:20122833. https://doi.org/10.1098/ 2013;280:1604–16. https://doi.org/10.1111/febs.12159. rspb.2012.2833. 29. Hundertmark M, Hincha DK. LEA (late embryogenesis abundant) proteins 9. Trujillo I, Ojeda MA, Urdiroz NM, Potter D, Barranco D, Rallo L, et al. and their encoding genes in Arabidopsis thaliana. BMC Genomics. 2008;9: Identification of the Worldwide Olive Germplasm Bank of Córdoba (Spain) 118. https://doi.org/10.1186/1471-2164-9-118. Julca et al. BMC Biology (2020) 18:148 Page 24 of 25 30. Durner J, Shah J, Klessig DF. Salicylic acid and disease resistance in plants. L.) in the Mediterranean Basin. Heredity. 2004;92:343–51. https://doi.org/10. Trends Plant Sci. 1997;2:266–74. https://doi.org/10.1016/S1360- 1038/sj.hdy.6800430. 1385(97)86349-2. 52. Doebley JF, Gaut BS, Smith BD. The molecular genetics of crop 31. Pearce G, Moura DS, Stratmann J, Ryan CA. RALF, a 5-kDa ubiquitous domestication. Cell. 2006;127:1309–21. https://doi.org/10.1016/j.cell.2006.12. polypeptide in plants, arrests root growth and development. Proc Natl Acad 006. Sci U S A. 2001;98:12843–7. https://doi.org/10.1073/pnas.201416998. 53. McKey D, Elias M, Pujol B, Duputié A. The evolutionary ecology of clonally 32. Vanholme R, Cesarino I, Rataj K, Xiao Y, Sundin L, Goeminne G, et al. propagated domesticated plants. New Phytol. 2010;186:318–32. https://doi. Caffeoyl shikimate esterase (CSE) is an enzyme in the lignin biosynthetic org/10.1111/j.1469-8137.2010.03210.x. pathway in Arabidopsis. Science. 2013;341:1103–6. https://doi.org/10.1126/ 54. Gaut BS, Díez CM, Morrell PL. Genomics and the contrasting dynamics of science.1241602. annual and perennial domestication. Trends Genet. 2015;31:709–19. https:// 33. Liu Q, Luo L, Zheng L. Lignins: biosynthesis and biological functions in doi.org/10.1016/j.tig.2015.10.002. plants. Int J Mol Sci. 2018;19 https://doi.org/10.3390/ijms19020335. 55. Terhorst J, Kamm JA, Song YS. Robust and scalable inference of population 34. Ashikawa I, Abe F, Nakamura S. DOG1-like genes in cereals: investigation of history from hundreds of unphased whole genomes. Nat Genet. 2017;49: their function by means of ectopic expression in Arabidopsis. Plant Sci. 303–9. https://doi.org/10.1038/ng.3748. 2013;208:1–9. https://doi.org/10.1016/j.plantsci.2013.03.011. 56. Klepo T, Toumi A, De La Rosa R, LeÓn L, Belaj A. Agronomic evaluation of 35. Rashid A, Deyholos MK. PELPK1 (At5g09530) contains a unique seedlings from crosses between the main Spanish olive cultivar ‘Picual’ and pentapeptide repeat and is a positive regulator of germination in two wild olive trees. J Hortic Sci Biotechnol. 2014;89:508–12. https://doi.org/ Arabidopsis thaliana. Plant Cell Rep. 2011;30:1735–45. https://doi.org/10. 10.1080/14620316.2014.11513113. 1007/s00299-011-1081-3. 57. Chu Y, Su X, Huang Q, Zhang X. Patterns of DNA sequence variation at 36. He S, Tan G, Liu Q, Huang K, Ren J, Zhang X, et al. The LSD1-interacting candidate gene loci in black poplar (Populus nigra L.) as revealed by single protein GILP is a LITAF domain protein that negatively regulates nucleotide polymorphisms. Genetica. 2009;137:141–50. https://doi.org/10. hypersensitive cell death in Arabidopsis. PLoS One. 2011;6:e18750. https:// 1007/s10709-009-9371-1. doi.org/10.1371/journal.pone.0018750. 58. Tuskan GA, Difazio S, Jansson S, Bohlmann J, Grigoriev I, Hellsten U, et al. 37. Lin R, Wang H. Arabidopsis FHY3/FAR1 gene family and distinct roles of its The genome of black cottonwood, Populus trichocarpa (Torr. & Gray). members in light control of Arabidopsis development. Plant Physiol. 2004; Science. 2006;313:1596–604. https://doi.org/10.1126/science.1128691. 136:4010–22. https://doi.org/10.1104/pp.104.052191. 59. Biswas S, Akey JM. Genomic insights into positive selection. Trends Genet. 38. Wang H, Wang H. Multifaceted roles of FHY3 and FAR1 in light signaling 2006;22:437–46. https://doi.org/10.1016/j.tig.2006.06.005. and beyond. Trends Plant Sci. 2015;20:453–61. https://doi.org/10.1016/j. 60. Pavlidis P, Živkovic D, Stamatakis A, Alachiotis N. SweeD: likelihood-based tplants.2015.04.003. detection of selective sweeps in thousands of genomes. Mol Biol Evol. 2013; 39. Mariotti R, Cultrera NGM, Díez CM, Baldoni L, Rubini A. Identification of new 30:2224–34. https://doi.org/10.1093/molbev/mst112. polymorphic regions and differentiation of cultivated olives (Olea europaea 61. Kaya HB, Akdemir D, Lozano R, Cetin O, Sozer Kaya H, Sahin M, et al. L.) through plastome sequence comparison. BMC Plant Biol. 2010;10:211. Genome wide association study of 5 agronomic traits in olive (Olea https://doi.org/10.1186/1471-2229-10-211. europaea L.). Sci Rep. 2019;9:18764. https://doi.org/10.1038/s41598-019- 40. Van de Paer C, Bouchez O, Besnard G. Prospects on the evolutionary 55338-w. mitogenomics of plants: a case study on the olive family (Oleaceae). Mol 62. Luu K, Bazin E, Blum MGB. pcadapt: an R package to perform genome scans Ecol Resour. 2018;18:407–23. https://doi.org/10.1111/1755-0998.12742. for selection based on principal component analysis. Mol Ecol Resour. 2017; 41. Huang J-L, Sun G-L, Zhang D-M. Molecular evolution and phylogeny of the 17:67–77. https://doi.org/10.1111/1755-0998.12592. angiosperm ycf2 gene. J Syst Evol. 2010;48:240–8. https://doi.org/10.1111/j. 63. Foll M, Gaggiotti O. A genome-scan method to identify selected loci appropriate 1759-6831.2010.00080.x. for both dominant and codominant markers: a Bayesian perspective. Genetics. 42. Drescher A, Ruf S, Calsa T, Carrer H, Bock R. The two largest chloroplast 2008;180:977–93. https://doi.org/10.1534/genetics.108.092221. genome-encoded open reading frames of higher plants are essential genes. 64. Burgarella C, Barnaud A, Kane NA, Jankowski F, Scarcelli N, Billot C, et al. Plant J. 2000;22:97–104. https://doi.org/10.1046/j.1365-313x.2000.00722.x. Adaptive introgression: an untapped evolutionary mechanism for crop 43. Wicke S, Schneeweiss GM, dePamphilis CW, Müller KF, Quandt D. The adaptation. Front Plant Sci. 2019;10:4. https://doi.org/10.3389/fpls.2019. evolution of the plastid chromosome in land plants: gene content, gene 00004. order, gene function. Plant Mol Biol. 2011;76:273–97. https://doi.org/10. 65. Zhou Y, Minio A, Massonnet M, Solares E, Lv Y, Beridze T, et al. The 1007/s11103-011-9762-4. population genetics of structural variants in grapevine domestication. Nat 44. Besnard G, Henry P, Wille L, Cooke D, Chapuis E. On the origin of the Plants. 2019;5:965–79. https://doi.org/10.1038/s41477-019-0507-8. invasive olives (Olea europaea L., Oleaceae). Heredity. 2007;99:608–19. 66. Baldoni L, Tosti N, Ricciolini C, Belaj A, Arcioni S, Pannelli G, et al. Genetic https://doi.org/10.1038/sj.hdy.6801037. structure of wild and cultivated olives in the Central Mediterranean basin. 45. Besnard G, Khadari B, Villemur P, Bervillé A. Cytoplasmic male sterility in the Ann Bot. 2006;98:935–42. https://doi.org/10.1093/aob/mcl178. olive (Olea europaea L.). Theor Appl Genet. 2000;100:1018–24. https://doi. 67. Davis PH. Flora of Turkey. Edinburgh: University Press; 1978. org/10.1007/s001220051383. 68. Lumaret R, Ouazzani N. Plant genetics. Ancient wild olives in Mediterranean 46. Fehrer J, Gemeinholzer B, Chrtek J, Bräutigam S. Incongruent plastid and forests Nature. 2001;413:700. https://doi.org/10.1038/35099680. nuclear DNA phylogenies reveal ancient intergeneric hybridization in 69. Besnard G, Terral J-F, Cornille A. On the origins and domestication of the Pilosella hawkweeds (Hieracium, Cichorieae, Asteraceae). Mol Phylogenet olive: a review and perspectives. Ann Bot. 2018;121:385–403. https://doi.org/ Evol. 2007;42:347–61. https://doi.org/10.1016/j.ympev.2006.07.004. 10.1093/aob/mcx145. 47. Barber JC, Finch CC, Francisco-Ortega J, Santos-Guerra A, Jansen RK. 70. Vargas P, Talavera S. Olea. Flora Ibérica. Madrid: CSIC; 2012. Hybridization in Macaronesian Sideritis (Lamiaceae): evidence from 71. Gaut BS, Seymour DK, Liu Q, Zhou Y. Demography and its effects on incongruence of multiple independent nuclear and chloroplast sequence genomic variation in crop domestication. Nat Plants. 2018;4:512–20. https:// datasets. Taxon. 2007;56:74–88. doi.org/10.1038/s41477-018-0210-1. 48. Pelser PB, Kennedy AH, Tepe EJ, Shidler JB, Nordenstam B, Kadereit JW, et al. 72. Zhou Y, Massonnet M, Sanjak JS, Cantu D, Gaut BS. Evolutionary genomics Patterns and causes of incongruence between plastid and nuclear of grape (Vitis vinifera ssp. vinifera) domestication. Proc Natl Acad Sci U S A. Senecioneae (Asteraceae) phylogenies. Am J Bot. 2010;97:856–73. https:// 2017;114:11715–20. https://doi.org/10.1073/pnas.1709257114. doi.org/10.3732/ajb.0900287. 73. Bartolini G, Petruccelli R. Classification, origin, diffusion and history of the 49. Green RE, Krause J, Briggs AW, Maricic T, Stenzel U, Kircher M, et al. A draft olive. Rome: FAO; 2002. sequence of the Neandertal genome. Science. 2010;328:710–22. https://doi. 74. Massei G, Hartley SE. Disarmed by domestication? Induced responses to org/10.1126/science.1188021. browsing in wild and cultivated olive. Oecologia. 2000;122:225–31. https:// 50. Malinsky M. Dsuite - fast D-statistics and related admixture evidence from doi.org/10.1007/PL00008850. VCF files. BioRxiv. 2019; https://doi.org/10.1101/634477. 75. Tang H, Zhang X, Miao C, Zhang J, Ming R, Schnable JC, et al. ALLMAPS: 51. Lumaret R, Ouazzani N, Michaud H, Vivier G, Deguilloux MF, Di Giusto F. robust scaffold ordering based on multiple maps. Genome Biol. 2015;16:3. Allozyme variation of oleaster populations (wild olive tree) (Olea europaea https://doi.org/10.1186/s13059-014-0573-1. Julca et al. BMC Biology (2020) 18:148 Page 25 of 25 76. Li H, Durbin R. Fast and accurate short read alignment with burrows- 98. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing wheeler transform. Bioinformatics. 2009;25:1754–60. https://doi.org/10.1093/ genomic features. Bioinformatics. 2010;26:841–2. https://doi.org/10.1093/ bioinformatics/btp324. bioinformatics/btq033. 77. Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a 99. Xie Z, Wang L, Wang L, Wang Z, Lu Z, Tian D, et al. Mutation rate analysis universal tool for annotation, visualization and analysis in functional via parent-progeny sequencing of the perennial peach. I. a low rate in genomics research. Bioinformatics. 2005;21:3674–6. https://doi.org/10.1093/ woody perennials and a higher mutagenicity in hybrids. Proc Biol Sci. 2016; bioinformatics/bti610. 283 https://doi.org/10.1098/rspb.2016.1016. 78. Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, et al. 100. Al-Shahrour F, Díaz-Uriarte R, Dopazo J. FatiGO: a web tool for finding Gapped BLAST and PSI-BLAST: a new generation of protein database search significant associations of gene ontology terms with groups of genes. programs. Nucleic Acids Res. 1997;25:3389–402. https://doi.org/10.1093/nar/ Bioinformatics. 2004;20:578–80. https://doi.org/10.1093/bioinformatics/ 25.17.3389. btg455. 79. Jones P, Binns D, Chang H-Y, Fraser M, Li W, McAnulla C, et al. InterProScan 101. Abascal F, Corvelo A, Cruz F, Villanueva-Cañas JL, Vlasova A, Marcet-Houben 5: genome-scale protein function classification. Bioinformatics. 2014;30: M, et al. Extreme genomic erosion after recurrent demographic bottlenecks 1236–40. https://doi.org/10.1093/bioinformatics/btu031. in the highly endangered Iberian lynx. Genome Biol. 2016;17:251. https:// 80. Waskom M, Botvinnik O, O’Kane D, Hobson P, Lukauskas S, Gemperline DC, doi.org/10.1186/s13059-016-1090-1. et al. Mwaskom/Seaborn: V0.8.1 (September 2017). Zenodo. 2017. doi: 102. McDonald JH, Kreitman M. Adaptive protein evolution at the Adh locus in https://doi.org/10.5281/zenodo.883859. drosophila. Nature. 1991;351:652–4. https://doi.org/10.1038/351652a0. 81. Dierckxsens N, Mardulyn P, Smits G. NOVOPlasty: de novo assembly of 103. Nielsen R, Williamson S, Kim Y, Hubisz MJ, Clark AG, Bustamante C. Genomic organelle genomes from whole genome data. Nucleic Acids Res. 2017;45: scans for selective sweeps using SNP data. Genome Res. 2005;15:1566–75. e18. https://doi.org/10.1093/nar/gkw955. https://doi.org/10.1101/gr.4252305. 82. Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, et al. SPAdes: a new genome assembly algorithm and its applications to single- cell sequencing. J Comput Biol. 2012;19:455–77. https://doi.org/10.1089/ Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in cmb.2012.0021. published maps and institutional affiliations. 83. Boetzer M, Henkel CV, Jansen HJ, Butler D, Pirovano W. Scaffolding pre- assembled contigs using SSPACE. Bioinformatics. 2011;27:578–9. https://doi. org/10.1093/bioinformatics/btq683. 84. Boetzer M, Pirovano W. Toward almost closed genomes with GapFiller. Genome Biol. 2012;13:R56. https://doi.org/10.1186/gb-2012-13-6-r56. 85. Slater GSC, Birney E. Automated generation of heuristics for biological sequence comparison. BMC Bioinformatics. 2005;6:31. https://doi.org/10. 1186/1471-2105-6-31. 86. Lowe TM, Chan PP. tRNAscan-SE on-line: integrating search and context for analysis of transfer RNA genes. Nucleic Acids Res. 2016;44:W54–7. https:// doi.org/10.1093/nar/gkw413. 87. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The genome analysis toolkit: a MapReduce framework for analyzing next- generation DNA sequencing data. Genome Res. 2010;20:1297–303. https:// doi.org/10.1101/gr.107524.110. 88. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27:2156–8. https://doi. org/10.1093/bioinformatics/btr330. 89. Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14:2611–20. https://doi.org/10.1111/j.1365-294X.2005.02553.x. 90. Earl DA, vonHoldt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012;4:359–61. https://doi.org/10.1007/s12686-011- 9548-7. 91. Rosenberg NA. DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes. 2004;4:137–8. https://doi.org/10.1046/j.1471-8286. 2003.00566.x. 92. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–75. https://doi.org/10.1086/ 519795. 93. Katoh K, Kuma K, Toh H, Miyata T. MAFFT version 5: improvement in accuracy of multiple sequence alignment. Nucleic Acids Res. 2005;33:511–8. https://doi.org/10.1093/nar/gki198. 94. Capella-Gutiérrez S, Silla-Martínez JM, Gabaldón T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics. 2009; 25:1972–3. https://doi.org/10.1093/bioinformatics/btp348. 95. Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post- analysis of large phylogenies. Bioinformatics. 2014;30:1312–3. https://doi.org/ 10.1093/bioinformatics/btu033. 96. Huson DH, Bryant D. Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. 2006;23:254–67. https://doi.org/10.1093/molbev/ msj030. 97. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995;57:289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x.