Genome sequence of the olive tree, Olea europaea

Background The Mediterranean olive tree (Olea europaea subsp. europaea) was one of the first trees to be domesticated and is currently of major agricultural importance in the Mediterranean region as the source of olive oil. The molecular bases underlying the phenotypic differences among domesticated cultivars, or between domesticated olive trees and their wild relatives, remain poorly understood. Both wild and cultivated olive trees have 46 chromosomes (2n). Findings A total of 543 Gb of raw DNA sequence from whole genome shotgun sequencing, and a fosmid library containing 155,000 clones from a 1,000+ year-old olive tree (cv. Farga) were generated by Illumina sequencing using different combinations of mate-pair and pair-end libraries. Assembly gave a final genome with a scaffold N50 of 443 kb, and a total length of 1.31 Gb, which represents 95 % of the estimated genome length (1.38 Gb). In addition, the associated fungus Aureobasidiumpullulans was partially sequenced. Genome annotation, assisted by RNA sequencing from leaf, root, and fruit tissues at various stages, resulted in 56,349 unique protein coding genes, suggesting recent genomic expansion. Genome completeness, as estimated using the CEGMA pipeline, reached 98.79 %. Conclusions The assembled draft genome of O. europaea will provide a valuable resource for the study of the evolution and domestication processes of this important tree, and allow determination of the genetic bases of key phenotypic traits. Moreover, it will enhance breeding programs and the formation of new varieties.

Spain and there are no legal restrictions on its use for research, including cv. Farga.
The tree age was estimated to be 1,200 years old based on dendrometric analyses (Antonio Prieto-Rodríguez personal communication). A combination of fosmid and whole genome shotgun (WGS) libraries were sequenced using Illumina sequencing equipment.
Primary data analysis was carried out using the standard Illumina pipeline (HCS 2.0.12.0, RTA 1.17.21.3). Mate-pair (MP) libraries (3,5,7 and 10 kb fragment sizes) were constructed at the CRG sequencing unit according to the Nextera Mate Pair Preparation protocol (Illumina Inc.), and sequenced on the HiSeq2500 platform in 2x150bp read length runs. The number of lanes and raw sequenced outputs for each library are summarized in Table 1.
Preliminary kmer analysis of PE data (Fig. 1) indicated a high level of heterozygosity in this sample. To reduce the risk of separately assembling two different haplotypes from the same locus and including them in the final assembly, a fosmid pooling strategy was chosen similar to the one used for the oyster genome project [2]. A fosmid library of 155,000 clones was constructed in the pNGS vector (Lucigen Corp.). Ninety-six pools of 1,600 clones each were made, and the purified DNA was used to construct short-insert PE libraries using the TruSeq™ DNA Sample Preparation Kit v2 (Illumina Inc.) and the KAPA Library Preparation kit (Kapa Biosystems) according to manufacturers' instructions. The pools were sequenced using TruSeq SBS Kit v3-HS (Illumina Inc.), in PE mode, 2×150 bp, in a fraction of a sequencing lane of the HiSeq2000 flowcell v3 (Illumina Inc.) according to standard Illumina operation procedures. The raw sequence yield per pool was 11.3 Gb on average (SD: 2 Gb), corresponding to~150 × depth. In addition a fosmid-end library was created from the same set of clones using the Lucigen pNGS protocol and run in one lane of a HiSeq2000. RNA was prepared from seven different tissues or developmental stages (root, young leaf, mature leaf, flower, flower bud, immature fruit, and green olives), using the Zymo ZR Plant RNA extraction kit (Zymo Research, Irvine, CA). Then, RNA-Seq libraries were prepared using the TruSeq™ RNA Sample Prep Kit v2 (Illumina Inc.) with minor modifications, and libraries were sequenced using the TruSeq SBS Kit v3-HS in PE mode with a read length of 2×75 bp. Over 50 million PE reads per sample were generated in a fraction of a sequencing lane on a HiSeq2000 (Illumina Inc.), following the manufacturer's protocol. Image analysis, base calling and run quality scoring were processed using the manufacturer's software Real Time Analysis (RTA 1.13.48), followed by generation of FASTQ sequence files using CASAVA software (Illumina Inc.).

Genome assembly
A kmer analysis was performed to estimate the genome size, level of heterozygosity and repeat content of the sequenced genome. Using the software Jellyfish v1.1.10 [3], 17-mers were extracted from the WGS PE reads (PE400), and unique kmers were counted and plotted according to kmer depth (Fig. 1). The homozygous or main peak is found at a depth of~52x. The estimated genome size (found by dividing the total number of kmers by the kmer depth of the main peak) is 1.38 Gb, which is at the low end of the range of empirical estimates. The C-value ranges from 1.45-2.33 pg (1.42 Gb-2.28 Gb), with the median at 1.59 pg (1.56 Gb) (data from [4], see [5][6][7][8][9]), suggesting the existence of variation in the repetitive fraction of the genome for the species. The left peak at 26x kmer depth indicates many polymorphic sites in the genome. In fact, using the Genomic Character Estimator program, gce v 1.0.0 [10], the heterozygous ratio based in kmer individuals is 0.054, and the corrected estimate of genome size is 1.32 Gb. Hereon the gce estimate is referred to as the 'assemblable' portion of the genome.
A pilot WGS assembly using only PE data was performed in order to generate enough contiguous sequences to gather library insert size statistics. PE reads were first filtered for contaminating sequences (phiX, Escherischia coli and other vector sequences, as well as O. europaea plastids) using GEM [11] with -m 0.02 (2 % mismatches). Then, the reads were assembled into scaffolds using AbySS v1.3.6 [12] with parameters: −s 600 − S 600-3000 − n 6 − N 10 − k 127 − l 75 − aligner map − q 10. This resulted in an assembly with a total length of 1.94 Gb, and contig and scaffold N50s of 3.7 kb and 3.8 kb, respectively. Library insert sizes were estimated by mapping against this draft assembly. For the WGS PE libraries sequenced on Illumina HiSeq2000 using 2x151 bp reads, the insert size distribution followed a bimodal distribution with a main peak at 725 bp and a smaller peak at 300 bp. Before continuing with the assembly, read pairs belonging to the smaller peak were filtered out, if connecting reads were found overlapping both mates of the pair. The inflated length (47 % of the assemblable part of the genome) and the poor contiguity obtained for the draft assembly are symptomatic of the expected difficulty in distinguishing divergent alleles of the same locus from true repeats. To address this challenge, the 96 sequenced fosmid pools (3.9x physical coverage of the genome, each pool covering~4 % of the genome) were assembled using the assembly pipeline shown in Fig. 2 to obtain 96 largely haploid assemblies (simulations of 1,600-clone pools with a genome size of 1.38 Gb show a mean of 2.5 % of sequenced bases to derive from separate overlapping clones, half of which would come from different alleles). Optimal kmer size was 97 for most of the pools. For each pool a base assembly was produced using ABySSv1.3.7 and parameters: −s 300 − S 300-5000 − n 9 − N 15 − k 97 − l 75 − aligner map − q 10. Afterwards, the base assemblies went through several rounds of gapfilling [13], decontamination, consistency checks, and rescaffolding with ABySSv1.3.7. The decontamination  Fig. 3), suggesting that a large fraction of fosmid clones are fully assembled. Mapping of fosmid ends to the merged assembly ('FP assembly' , see below) gives an estimate of the clone insert size distribution (mean of 36.7 kb ± SD 4.97 kb) that corresponds well with the right peak of the scaffold sizes. The 96 fosmid pool assemblies were then merged based on overlaps using in-house OLC-like assemblymerging software called ASM (L. Frias and P. Ribeca, manuscript in preparation; scripts are publicly available at [14]. Two rounds of merging were performed, with intermediate scaffolding and gapfilling steps. In the first round, a minimum overlap of 2,400 bp and high sequence similarity (maximum edit distance of 1.5 %) was used, while in the second round, longer overlaps (4,000 bp) and higher sequence divergence (maximum edit distance of 10 %) were used in order to merge allelic regions. Each round of merging collapses repeats unless higher order information supports a unique path for resolving a repetitive region; this includes both the sequence of the input data (contigs) and scaffolding information (i.e., the order of contigs in scaffolds in the original fosmid pool assemblies). Merging produced an intermediate assembly (named 'FP assembly' in Fig. 4) with a scaffold N50 of~45 kb and a total length of 1.38 Gb. Although this assembly was 4.54 % larger than the assemblable genome size (1.32 Gb), gene completeness according to CEGMA was only 95.97 % complete and 97.58 % partial, suggesting that 2.42-4.03 % of the gene space may have been missed.
To increase the overall completeness of the assembly, all WGS reads that did not map to the FP assembly were selected and used to obtain a complementary assembly using ABySSv1.5.2 with parameters: −s 300 − S 300-5000 − n 10 − N 10 − k 95 − l 75 − aligner map − q 10. This assembly accounts for 60.7 Mbp of sequence, and has an N50 of 1,506 bp for contigs and 2,351 bp for scaffolds. This assembly was then broken into contigs, 50 bp was eroded from the ends of each contig, then contigs smaller than 200 bp were filtered out. Both assemblies were subsequently gathered by joining the WGS contigs with the merged fosmid pool assembly, and scaffolding them with SSPACE 2.0 [15]. To account for read pairs coming from two different alleles in the same genomic region, reads were mapped to the SSPACE input assembly with gem-mapper (settings: m = 0.05 and e = 0.1) and filters were applied to detect unique mappings with no subdominant match. The resulting comprehensive assembly had a scaffold N50 of 303.7 kb and a total length . The excess of assembled sequence is likely to be caused by the presence of artificial duplications during the assembly process (i.e., uncollapsed haplotypes that have been resolved in two different contigs). Several strategies were used to refine the assembly and obtain a haploid reference. First, consistency check was applied to remove local misassemblies by mapping short and intermediate libraries (PE720, MP3k and MP5k) to the input assembly: a positive score is assigned to the assembly regions supported by read pairs separated by distances falling within the limits (mean ± 3σ) of the empirical distribution, while a negative score is assigned to regions where read pairs map i) outside of these bounds, ii) in inconsistent orientation, or iii) to different scaffolds. Regions where the sum of these two vectors is negative are removed from the assembly. After applying this consistency check, the resulting assembly had 46,893 consistent contig blocks (compared to 25,042 contigs before the consistency check), giving a total of 1.46 Gb with an N50 of 101 kb. Second, this assembly was collapsed using a minimum overlap of 4 kb and the gem-mapper parameters − e 0.03 and − m 0.02, so only close matches were merged (similar uncollapsed haplotypes, identical assembly artifacts, and near identical repeats). Additionally, in order to avoid spurious joins, tip merging was applied to the alignment graph down to overlaps of 250 bp. Finally, no repeat resolution was applied, but coherent links from input scaffolds were reinserted. Consequently, the assembly length shrunk to~1.30 Gb, almost matching the assemblable fraction of the genome (1.32 Gb). An additional consistency check was run on the collapsed assembly using the short and intermediate libraries (PE720, MP3k and MP5k), which resulted in breaking the assembly from 64,814 into 72,593 scaffolds, giving a total length of 1.30 Gb with a scaffold N50 of 50 kb. This assembly length is what was expected based on the gce estimate. As a final assembly step, PE reads with high divergence (gem-mapper parameters m = 0.05 and e = 0.08) were mapped to the assembly and rescaffolded with SSPACE 2.0 using parameters k = 3 and a = 0.6. Then, scaffolds shorter than 500 bp were discarded, and the GapFiller program [13] was used to close about 40 % of the assembly gaps. This assembly was labeled 'Oe3'.
The Oe3 assembly was polished using a mappingbased strategy designed to correct single nucleotide substitution and short insertion-deletion errors. First, one library of paired-end reads (PE725) was aligned using BWA mem (v0.7.7) [16] and variant calling was performed. Selecting only homozygous alternative variants, an alternative FASTA sequence was obtained using GATK (v3.5) FastaAlternateReferenceMaker [17]. After discarding scaffolds shorter than 500 bp, the resulting assembly (Oe5) had a scaffold N50 of 444 kb and a contig N50 of 51 kb. After detecting putative contamination in some scaffolds of the Oe5 assembly, a final decontamination step was performed against yeast, bacteria, arthropod and mitochondrial sequences, combining homology search results obtained by BLAST and, in the case of mitochondrial sequences, regions of high depth (~6000x). In total, 509 scaffolds were deleted from Oe5 and some parts of another 27 scaffolds were removed. The assembly resulting from this step, Oe6, has a scaffold N50 of 443 kb and a contig N50 of 52 kb (Table 2). Oe6 contains 48,419 gaps comprising 53,969,601 sites. The gene completeness of this assembly was estimated using CEGMA [18] and BUSCO (Benchmarking Universal Single-Copy Orthologs) [19]. CEGMA analysis resulted in a gene completeness of 98.79 %, while BUSCO, using a plant-specific database of 956 genes, determined a completeness of 95.6 % of plant genes. A summary of the complete assembly strategy is shown in Fig. 4.
Partial assembly of an olive tree associated fungus: Aureobasisium pullulans One of the putative sources of non-plant sequence present in the olive samples was considered of interest; it was represented among the fosmid pools and seemed to belong to the fungal genus Aureobasidium, which has been previously associated with olive trees [20]. To assemble a partial sequence of this genome, four fully sequenced Aureobasidium genomes were downloaded from JGI [21]. Then, BWA v0.7.3a [16] was used to map all the reads from the fosmid libraries to the four genomes. Once mapped, the reads were filtered allowing only soft clipping for a maximum of one-third of the read, and deleting read pairs when only one of the pairs passed the filters. This resulted in a collection of 18,549,090 reads, which were assembled with SPAdes (See figure on previous page.) Fig. 3 Fosmid pool assembly pipeline. For each fosmid pool, a single paired-end (PE) library sequenced at 2 x 150 bp was first filtered and trimmed of pNGS vector sequences, as well as those of Escherichia coli and other common contaminants, including Olea europaea chloroplast sequences. Reads were assembled with ABySS, gapfilled with GapFiller, and contaminants removed using a BLAST homology search. A consistency check was performed, breaking the assemblies at any point inconsistent with the proper insert size and orientation of fosmid pool PE reads. The resulting contigs were scaffolded using whole genome shotgun (WGS) data, followed by another round of gapfilling, decontamination and consistency checking, this time including the new WGS data. To repair the consistency broken assembly, a final round of scaffolding, gapfilling and decontamination was performed v.3.1.1 [22]. Scaffolding was done using the assembled fosmids using SSPACE-LongRead [23], and gaps were filled with gapcloser [24]. These two steps were repeated twice. The final alignment was then compared to the Aureobasidium genomes using BLAST. Contigs longer than 200 nt, for which less than 20 % of their sequence mapped against any of the Aureobasidium genomes, were separated and compared against the NCBI nonredundant nucleotide database [25]. Only those contigs with first hits to fungal species were kept. The final assembly comprised 18 Mb, roughly two-thirds of the typical size of Aureobasidium genomes (25)(26)(27)(28)(29). To identify the species and strain, the most common fungal markers used for fungal barcoding were identified (ITS, SSU, LSU, RPB1, RPB2 and EF1). Most of the markers were missing in the assembly or were too short; based on a 769 nt fragment of the RPB1 gene, the most similar sequence was that of Aureobasidium pullulans isolate AFTOL-ID 912 (DQ471148.1); a strain that was isolated from the grape plant Vitis vinifera. The identity of this fragment was 99.95 % indicating that this was likely a different strain of the same species. Augustus [26] was used to perform gene annotation. The training parameters were obtained using scaffold 1 of the published A. pullulans genome, and then used to predict proteins in our strain of A. pullulans. This resulted in 6,411 proteins.

Olive tree genome annotation
To annotate the olive tree genome, consensus gene models were obtained by combining transcript alignments, protein alignments, and gene predictions. A flowchart outlining these steps is shown in Fig. 5. Transcripts for assembly with Program to Assemble Spliced Alignments (PASA; r2014-04-17) [27] were obtained as follows: first, RNA-Seq reads generated from different tissues by our group (see above), plus publicly available datasets in the Sequence Read Archive (SRA) ( Table 3), were aligned to the final assembly Oe6 with GEM v1.6.1 [11]. Transcript models were subsequently generated using the standard Cufflinks v2.1.1 pipeline [28] starting with the BAM files, resulting in 2,056,606 transcripts, which were then added to the PASA database. In addition, 12,959 olive expressed sequence tags (ESTs) and mRNAs present in Genbank (October 27, 2014) [29][30][31] were also added to PASA using GMAP v2013-10-28 [32] as the alignment engine. All of the above transcript alignments were then assembled by PASA, resulting in 942,302 PASA assembled transcripts, which were scanned with PASA's Transdecoder program [27] to detect likely protein coding regions. This tool predicted a total of 169,562 candidate genes. From these, a training set for ab initio gene predictors was created from PASA models coding for complete proteins, longer than 500 amino acids and with a BLAST hit to either the Lamiidae or Asteridae proteomes. A training set of 589 non-  For ab initio gene prediction, transposable element repeats in the Oe6 assembly were first masked with RepeatMasker v4-0-5 [34] using a custom repeat library constructed by running RepeatModeler v1-0-7 and adding some olive-specific repeats [35]. A search was also carried out for masked proteins encoded by transposable elements (TEs) provided in the RepeatMasker Library of TE proteins. Low complexity repeats were left unmasked for this purpose. In total, 63 % of the assembly was masked.
On this masked assembly four different ab initio gene predictors were run, since combiners like EvidenceModeler work better when finding consensus among the output of a diverse set of gene prediction algorithms, and orthogonal evidence such as transcript and protein mapping. O. europaea protein-coding gene predictions were obtained with GeneID v1.4.4 [36] trained specifically for O. europaea with GeneidTrainer using the training set of 589 genes; with Augustus v3.0.2 [26] trained with the etraining script that comes with Augustus using the same training set; and with GlimmerHMM v3.0.1 [37] trained with the trainGlimmerHMM script that comes with the program using the same training set. Finally, GeneMark-ES v2.3 [38] gene predictions were obtained by running it in its self-trained mode. The number of predicted gene models ranged from 48,237 with GeneMark-ES to 97,542 with GlimmerHMM. Geneid, Augustus and Genemark-ET v4.21 were also used to generate predictions incorporating intron evidence, which was extracted from the RNA-Seq data, by obtaining the junctions after mapping it with GEM (see below). Junctions overlapping with ab initio GeneID predictions, Augustus predictions, or with protein mappings were taken as intron evidence. Running GeneID with hints resulted in a total set of 74,231 gene models; Augustus with hints resulted in 70,906; and Genemark-ET with 64,329 gene models.
Evidence Modeler r2012-06-25 (EVM) [39] was used to obtain consensus CDS models using the three main sources of evidence described above: gene predictions, aligned transcripts and aligned proteins. EVM was run with three different sets of evidence weights, and the resulting consensus models with the best specificity and sensitivity as determined by intersection (BEDTools v2.16.2 intersect [40]) with the transcript mappings, were chosen for the final annotation (Table 4 shows the best-performing weights). Consensus CDS models were then updated with untranslated regions (UTRs) and alternative exons through two rounds of PASA annotation updates. A final quality control was performed to fix reading frames and intron phases, and remove some transcripts predicted to be subject to nonsense-mediated decay. The resulting transcripts were clustered into genes using shared splice sites or substantial sequence overlap as criteria for designation as the same gene. This resulted in a preliminary set of 56,349 proteincoding genes, whose 89,982 transcripts encode 79,910 unique protein products (~1.59 transcripts per gene). Systematic identifiers with the prefix 'OE6A' were assigned to the genes, transcripts and derived protein products. Functional annotation was performed with InterProScan-5.17-56.0 [41], 30,900 protein-coding genes were annotated with gene ontology (GO) terms, and 41,257 were assigned a function.
The predicted O. europaea protein-coding set was then compared with those in four other selected plant genomes (Arabidopsis thaliana, Erythranthe guttata, Solanum lycopersicum, and Ricinus communis) downloaded from the   [42] with an e-value less than 0.01 and with at least 50 % of identity (Table 5). General statistics for transcript, coding sequence and exon lengths in O. europaea are similar to those in the other species, but the number of genes is significantly larger. The number of exons per transcript is slightly lower than in the four compared species. It is possible that more false-positive single-exon genes have been annotated; however, the number of single-exon CDS is not higher, although there is a slight shift in the distribution toward fewer coding exons per transcript (Fig. 6).
The increased number of coding genes in O. europaea suggests the existence of a large-scale genome duplication with respect to the other species. Although this possibility deserves more detailed analysis, preliminary analyses of gene comparisons identified 34,195 O. europaea genes with O. europaea paralogs that are more similar to each other than to the corresponding best hit in E. guttata (80.5 % of the total proteins with hits in E. guttata), the closest species in this analyses. Also, from the 14,437 paralogous pairs found in O. europaea that represent each other's reciprocal best hit, 10,711 pairs had the same best hit in E. guttata (which represents 74.2 % of the pairs). These results suggest that a high  proportion of the O. europaea gene repertoire has been duplicated since the separation of these two lamiales species. To discard the possibility that these duplicates resulted from uncollapsed heterozygous alleles, heterozygous single nucleotide variants (SNVs) identified by variant calling using samtools mpileup in pairs of putatively recent duplicates were counted and compared with those in singletons (genes without recent paralogs). The mean is significantly higher in genes within recent duplicate pairs (Welch's Two Sample t-test p-value < 2.2e-16). Finally, the 70 % quantile of two-copy SNV counts is 42 and 8 for the one-copy genes. In the case where uncollapsed (duplicated) alleles are frequent, one would expect to obtain the opposite pattern, as reads coming from the same locus would independently map to one of the two uncollapsed haplotypes in the assembly, thus dramatically reducing the number of heterozygous SNVs called. Although further and more detailed analyses are required, these results suggest extensive gene duplication in the lineage leading to the olive tree. The possibility of a whole genome duplication is consistent with the increased chromosomal number in O. europaea (2n = 46), as compared to closely related lamiales such as Erythranthe guttata (2n = 28) [43] and Sesamum indicum (2n = 26) [44]. Non-coding RNAs (ncRNAs) were annotated by running the following steps. First, the program cmsearch (v1.1) that comes with Infernal [45] was run with the Rfam database of RNA families (v12.0) [46]. Also, tRNAscan-SE (v1.23) [47] was run in order to detect the transfer RNA genes present in the genome assembly. To detect long non-coding RNAs (lncRNAs), PASA assemblies that had not been included in the annotation of protein-coding genes (i.e., expressed genes that were not translated to protein) were first selected. Those longer than 200 bp and with a length not covered by a small ncRNA at least 80 % were incorporated into the ncRNA annotation as lncRNAs. The resulting transcripts were clustered into genes using shared splice sites or significant sequence overlap as criteria for designation as the same gene. Systematic identifiers with the prefix 'OE6ncA' were assigned to the genes and their derived transcripts. In total, 25,199 non-coding genes have been annotated, among which 20,082 are lncRNAs.
In summary, we report the first genome sequencing, assembly, and annotation of the Mediterranean olive tree. This genome assembly will provide a valuable resource for studying developmental and physiological processes, investigating the past history of domestication, and improving the molecular breeding of this economically important tree.