Gautier ; Sarah

Félix Homa

Stéphanie ; Pointet

Sandy ; Contreras

François ; Sabot

Benoit Nabholz

Sylvain S Santoni

Laure ; Sauné

Morgane ; Ardisson

Nathalie N Chantret

Felix ; Homa

Nabholz

; Benoit

Sylvain ; Santoni

Nathalie ; Chantret

Christopher ; Sauvage

James ; Tregear

Cyril ; Jourda

David ; Pot

Yves Vigouroux

Felix Homa

Sylvain Santoni

Nathalie Chantret

Hana Chair

Nora Scarcelli

Claire Billot

Nabila Yahiaoui

Roberto Bacilieri

Bouchaib Khadari

Michel Boccara

Adéline Barnaud

Jean-Pierre Péros

Jean-Pierre Labouisse

Jean-Louis Pham

Jacques David

Sylvain Glémin

Manuel Ruiz
email: manuel.ruiz@cirad.fr

Manuel Ruiz

A large set of 26 new reference transcriptomes dedicated to comparative population genomics in crops and wild relatives

Keywords: Reference transcriptomes, RNA-Seq, Comparative genomics

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Population genomic approaches including both crops and wild relatives are essential to provide key information on the crop domestication processes and the dynamics of recent adaptation [START_REF] Ross-Ibarra | Plant domestication, a unique opportunity to identify the genetic basis of adaptation[END_REF]. This has been applied to an increasing number of species but mainly annual crops, and especially cereals [START_REF] Hufford | Comparative population genomics of maize domestication and improvement[END_REF]. However, domestication scenarios may vary between species and more insights could be gained by comparative approaches in diverse taxa [START_REF] Glemin | A comparative view of the evolution of grasses under domestication[END_REF][START_REF] Meyer | Patterns and processes in crop domestication: an historical review and quantitative analysis of 203 global food crops[END_REF]. It could be also very fruitful to compare molecular evolutionary patterns among species with contrasted life-history or ecological traits. Life-history or ecological traits may influence genome evolution through their effect on key population genetic parameters (effective size, recombination rates, and mutation rates).

Furthermore, studying the evolution of gene families and relating it to expression data across lineages may help to identify which molecular functions play a key role in adaptation. A crucial initial step for such studies is to obtain a large set of reference genes in different crops, both for wild and cultivated populations, but also in close relative species as much insight can be gained by comparing genomes beyond the domestication scale; for instance several selection detection methods require outgroup species [START_REF] Vitti | Detecting natural selection in genomic data[END_REF]. Today, the massive sequencing technologies can produce essential resources for large comparative genomic studies among many species, including non-model species [START_REF] Ellegren | Genome sequencing and population genomics in non-model organisms[END_REF].

We have made available a unique dataset of 26 new plant reference transcriptomes that were produced, assembled and annotated with the same protocol. We chose a set of various Mediterranean and tropical crop species with contrasted life-history traits (ex: annual herbs/trees, outcrossers/selfers) and with a wide phylogenetic distribution (including both monocots and eudicots belonging to different orders). For non-model

Accepted Article

This article is protected by copyright. All rights reserved.

crops (ex: yam, fonio millet), we chose both the crop and close relatives. For some species already having a reference genome (ex: tomato, cocoa) we focused on relative outgroup species (see detailed list in Table 1). We chose high-throughput RNA-Seq technology for feasibility reasons, for possibilities of obtaining the coding part of the genome and expression levels, and for comparing the fate of some genes in different species. As we were focusing on providing tools for comparative and population genomics analyses, we were particularly attentive to the optimization of the protocol for assembling the Illumina short-reads in order to reconstruct the best ratio of unambiguous correct cDNAs, which we favoured instead of the exhaustive representation of transcriptomes.

We assessed and compared the quality of de novo assemblies using data from two crops for which an annotated reference genome is available, grapevine and sorghum.

This study allowed us to provide a useful set of plant transcriptomes, contributing both to the enrichment of genomic resources for various crops of agronomic interest (some of them having being neglected so far such as yam, fonio millet or einkorn wheat) and phylogenetically related species that can be used as outgroups, which are invaluable tools for many evolutionary studies (Table 1). Among the chosen crops, some of them have been neglected so far while they can contribute to taking up some current agronomic challenges. For example, the diploid einkorn wheat could provide a useful resource for polyploidy wheat breeding programs, underutilized African cereals such as fonio millet deserve to be developed as an alternative crop well adapted to tropical conditions, and yam is a good model to understand clonally propagated crops.

Accepted Article

This article is protected by copyright. All rights reserved.

Materials and Methods

Preparation of RNA samples

Samples were constituted by combining different organs, including leaves, fruits/grains, and inflorescence tissues. Samples were ground in liquid nitrogen and total cellular RNA was extracted using a Spectrum Plant Total RNA kit (Sigma, Inc., USA) with a DNAse treatment. Oil palm total RNAs were extracted in a similar way using an RNeasy kit (Qiagen) and Mauritia flexuosa total RNAs were extracted as described by [START_REF] Morcillo | Somaclonal variation in micropropagated oil palm. Characterization of two novel genes with enhanced expression in epigenetically abnormal cell lines and in response to auxin[END_REF]. RNA concentration was first measured using a NanoDrop ND-1000 Spectrophotometer then with the Quant-iT™ RiboGreen® (invitrogen) protocol on a Tecan Genius spectrofluorimeter. RNA quality was assessed by running 1 μL of each RNA sample on RNA 6000 Pico chip on a Bioanalyzer 2100 (Agilent Technologies, Inc., USA). Samples with an RNA Integrity Number (RIN) value greater than eight were deemed acceptable according to the Illumina TruRNA-Seq protocol. Usually, for each genotype, a 2 μg mixture consisting of RNA from the inflorescence (80%) and leaves (20%) was prepared. For a few species, three sources of RNA were used: 65% from the inflorescence, 15 % from leaves, and 20% from fruits or grains when these additional organs were of specific agronomic interest and could increase the transcriptome representation.

Illumina library production

The TruSeq RNA sample Preparation v2 kit (Illumina Inc., USA) was used according to the manufacturer's protocol with the following modifications. Poly-A containing mRNA molecules were purified from 2 ug of total RNA using poly-T oligo attached magnetic beads. The purified mRNA was fragmented by addition of the fragmentation buffer and was heated at 94°C in a thermocycler for 4 min. A fragmentation time of 4 min was

Accepted Article

This article is protected by copyright. All rights reserved.

used to yield library fragments of 250-300 bp. First strand cDNA was synthesized using random primers to eliminate the general bias towards the 3' end of the transcript.

Second strand cDNA synthesis, end repair, A-tailing, and adapter ligation was performed in accordance with the protocols supplied by the manufacturer. Purified cDNA templates were enriched by 15 cycles of PCR for 10 s at 98°C, 30 s at 65°C, and 30 s at 72°C using PE1.0 and PE2.0 primers, and with Phusion DNA polymerase. Each indexed cDNA library was verified and quantified using a DNA 100 Chip on a Bioanalyzer 2100 then equally mixed by ten (from different genotypes). The final library was then quantified by real time PCR with the KAPA Library Quantification Kit for Illumina Sequencing Platforms (Kapa Biosystems Ltd, SA) and adjusted to 10 nM in water prior to sequencing at the Montpellier Genomix platform (http://www.mgx.cnrs.fr/).

Illumina library clustering and sequencing conditions

The final mixed cDNA library was sequenced using the Illumina mRNA-Seq, paired-end protocol on a HiSeq2000 sequencer, for 2 x 100 cycles. The library was diluted to 2 nM with NaOH. 2.5 μL was transferred into 497.5 μL HT1 to give a final concentration of 10 pM. 120 μL was then transferred to a 200 μL strip tube and placed on ice before loading onto the Cluster Station. A mixed library consisting of 10 individually indexed libraries was run on a single lane. Flow cells were clustered using Paired-End Cluster Generation Kit V4, following the Illumina PE amplification Linearization Blocking PrimerHyb v7 recipe. Then, the flow cell was loaded onto the Illumina HiSeq 2000 instrument following the manufacturer's instructions. The sequencing chemistry used was v4 (FC-104-4001, Illumina) using SCS 2.6 and RTA 1.6 software with the 2 x 100 cycles, paired-end, indexed protocol. Illumina base calling files were processed using

Accepted Article

This article is protected by copyright. All rights reserved.

the GERALD pipeline to produce paired sequence files containing reads for each sample in Illumina FASTQ format.

Read cleaning

Reads were preprocessed with cutadapt (Martin 2011) using the TruSeq index sequence corresponding to the sample, searching within the whole sequence. We also used cutadapt to trim the end of the reads with low quality scores (parameter -q 20) and kept reads with a minimum length of 35 bp. We then filtered the reads on the basis of their mean quality score, keeping those with a mean quality higher than 30.

Subsequently we separated orphan reads (i.e those for which the mate was discarded in the previous steps) using a homemade script.

Assembly protocols

Reads were assembled using Abyss [START_REF] Simpson | ABySS: a parallel assembler for short read sequence data[END_REF]) followed by one step of Cap3 [START_REF] Huang | CAP3: A DNA sequence assembly program[END_REF]. Reads returned as singletons by the first assembly run were discarded. For assembly, only paired reads were used. Abyss was launched using the paired-end option with a kmer value of 60. Cap3 was launched with the default parameters, notably 40 bases of overlap and the percentage of identity was set at 90%.

Comparison of assembly protocols

We compared different combinations using the assemblers Abyss [START_REF] Simpson | ABySS: a parallel assembler for short read sequence data[END_REF]), Cap3 [START_REF] Huang | CAP3: A DNA sequence assembly program[END_REF], and Trinity [START_REF] Haas | De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis[END_REF], with different parameters. We compared assemblies with Trinity (T), Abyss (A), and Abyss followed by one step of Cap3 (A/C), and Abyss followed by two steps of Cap3 (A/C/C) according to the strategy of [START_REF] Cahais | Reference-free transcriptome assembly in non-model animals from next-generation sequencing data[END_REF]. We used these protocols with reads from either only one individual or from multiple individuals. When the assemblies with only one

Accepted Article

This article is protected by copyright. All rights reserved.

individual were used, we chose the individual with the highest read number after cleaning: EC4 for Sorghum bicolor (SSM1123) and VC2 for Vitis vinifera (Morenoa).

When assemblies were performed with reads from multiple individuals, we used 10 individuals for Sorghum bicolor and 12 individuals for Vitis vinifera (Table 2). In both cases, all reads came from a single sequencing lane. Since sampling, bank preparation and sequencing conditions were the same for all individuals, we assumed that the individual with the highest read number would enable reconstruction of the maximum number of weakly expressed transcripts. To perform Trinity with reads from multiple individuals, we used in silico read normalization to reduce memory requirements and improve computational times. Trinity was launched with the default parameters, except that the threshold for minimum contig length was increased from 200 bp to 250 bp.

Criteria for comparing the quality of assemblies

For each assembly we calculated the number of contigs, the mean length of the contigs, the N50, and the number of contigs that mapped to the reference genome (see Exonerate mapping description below). We used BWA [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] (3 mismatches allowed) to map the reads from 10 individuals for Sorghum bicolor, and 12 individuals for Vitis vinifera to the corresponding assessed contigs. For each assembly, we calculated the percentage of mapped reads and the percentage of correctly paired mapped reads. Correctly paired mapped reads were defined as being within 1000 bases of each other. We also mapped the reads to the transcriptome predicted from the reference genome annotations.

Reference transcriptomes

For Vitis vinifera, the sequences of the transcriptome derived from the reference genome annotation were downloaded from the CRIBI genomes portal (http://genomes.cribi.unipd.it/DATA/V1/FASTA/V1_mrna.fa).

Accepted Article

This article is protected by copyright. All rights reserved.

The Sorghum bicolor reference transcriptome corresponds to the Sorghum bicolor v1.4 dataset available on the phytozome web site (http://www.phytozome.net/) at ftp://ftp.jgipsf.org/pub/compgen/phytozome/v9.0/Sbicolor_v1.4/annotation/ (Sbicolor_79_transcript.fa).

BLAST analysis of contigs against the corresponding transcriptome predicted from the genome

For each contig, a BLASTN search was performed against the appropriate transcriptome reference. We only considered the first ten hits with an e-value of 1e-5.

Hits were considered significant when (i) the alignment length (merging all high-scoring segment pairs) was at least 50% of the query sequence or at least 50% of the hit sequence, and (ii) sequence identity between query and hit was more than 90% across the alignment. This last threshold was chosen higher compared to [START_REF] Cahais | Reference-free transcriptome assembly in non-model animals from next-generation sequencing data[END_REF] because contigs were here compared to the reference transcriptome of the same species.

Classification of BLAST results

Contigs were classified into seven categories, defined by the number and nature of significant BLAST hits. Based on [START_REF] Cahais | Reference-free transcriptome assembly in non-model animals from next-generation sequencing data[END_REF] we distinguished the following cases: a) Contigs with no significant hit were called 'No Hit'. b) Contigs with a single significant hit and such that no other contig found the same hit. Within this category, contigs were called 'Full length' if the contig ⁄ hit alignment covered at least 90% of the hit sequence or 'Partial' otherwise. c) Contigs with a single significant hit shared by other contigs.

Accepted Article

This article is protected by copyright. All rights reserved.

Within this category, when two contigs overlapped, if the length of the hit sequence to which both contig queries aligned covered at least 50% of the shortest of the two contigs, contigs were called "Allele", whereas those showing no significant overlap with any other contigs were called "Fragment". d) Contigs with several significant hits, all specific to this contig.

When the hit sequences showed significant overlap, with the hit/hit alignment covering at least 50% of the shortest of the two hit sequences, contigs were called "Multi copies". Whereas contigs whose hits showed no significant overlap were called "Chimera". e) Contigs not included in any of the above categories, involving more complex patterns.

If m contigs had n significant hits, and if m and n were equal, and when the m first hits of the m contigs were distinct from each other, contigs were called 'Full length' or 'Fragment', depending on the contig ⁄ hit alignment. In all other cases, contigs were called 'Multi copies'.

Mapping of the contigs on the reference genome

We used Exonerate [START_REF] Slater | Automated generation of heuristics for biological sequence comparison[END_REF] to perform mapping of the contigs on the reference genome. We used the following parameters: sizes of introns limited to 5,000 bp, model set to est2genome, and percent 50 and bestn 5 to keep the best five hits if their scores were at least 50% of the best hit. In order to easily compare the positions of our mapping with the predicted mRNAs, we used the parameter showtargetgff which gave the alignment in GFF format.

For the percentage of annotations mapped by contigs, we took into account the number of distinct annotations mapped. The annotations corresponded to the

Accepted Article

This article is protected by copyright. All rights reserved.

determined positions of mRNAs on the genome. The number of mapped annotations was determined by the tool IntersectBed from bedtools [START_REF] Quinlan | BEDTools: a flexible suite of utilities for comparing genomic features[END_REF].

BBMH (Best Blast Mutual Hits) comparison of transcriptomes

For Solanum melongena, the sequences of the transcriptome derived from annotation of the reference genome SME_r2.5. Within each transcriptome, we used CD-HIT [START_REF] Fu | CD-HIT: accelerated for clustering the next-generation sequencing data[END_REF] to remove highly similar sequences. To access pairs of similar sequences between compared data sets, Best Blast Mutual Hits (BBMH) analyses were computed with BLASTN and an e-value of 1e-10. Venn diagrams were created with jvenn [START_REF] Bardou | jvenn: an interactive Venn diagram viewer[END_REF].

Functional Annotations

Contigs from the different assemblies were blasted as follows: firstly, they were compared to the Swissprot database. Queries with no hit were compared to the TrEMBL database. Then, queries with no hit against the TrEMBL database were compared to the NR database, and sequences still having no hit were compared against the NT database. We considered hits with an e-value of 1e -5 and only the best five hits were kept.

Sequences were translated using prot4EST [START_REF] Wasmuth | prot4EST: translating expressed sequence tags from neglected genomes[END_REF]). We provided it with the outputs of the proteic BLAST (against Swissprot, TrEMBL and NR).

Accepted Article

This article is protected by copyright. All rights reserved.

The BLAST output files were loaded into Blast2GO [START_REF] Gotz | High-throughput functional annotation and data mining with the Blast2GO suite[END_REF] to perform Gene Ontology (GO) annotation and GO terms enrichment analysis. We used the three GO ontologies: molecular function, biological process, and cellular component. We categorized the enriched GO terms with CateGOrizer (Hu et al. 2008) using the GO_slim2 classification method and single occurrences count.

Results and Discussion

Validation of the assembly protocol

The performance of the different computational methods available for transcript reconstruction from RNA-Seq data depends on various criteria, including: the type of organism to be assembled; the preference given to either sensitivity or precision of the predictions; and the type of analysis to be performed with the reconstructed contigs [START_REF] Steijger | Assessment of transcript reconstruction methods for RNA-seq[END_REF]. Here, the criterion was to maximize production of accurately predicted gene sets that can be commonly found in the transcriptomes of different individuals of a species for intra-species analyses (e.g. SNP detection, population genomics). We followed the overall strategy of [START_REF] Cahais | Reference-free transcriptome assembly in non-model animals from next-generation sequencing data[END_REF] because their protocol focused on optimizing the reconstruction of unambiguous correct cDNAs for comparative and population genomics purposes. We refined this assembly protocol by re-evaluating the utility of combining the assemblers Abyss [START_REF] Simpson | ABySS: a parallel assembler for short read sequence data[END_REF]) and Cap3 [START_REF] Huang | CAP3: A DNA sequence assembly program[END_REF]. We assessed and compared the quality of de novo Illumina short-read assemblies using data from two crops for which an annotated reference genome is available, grapevine (Vitis vinifera subsp. vinifera) [START_REF] Jaillon | The grapevine genome sequence suggests ancestral hexaploidization in major angiosperm phyla[END_REF]) and sorghum (Sorghum bicolor subsp. bicolor) [START_REF] Paterson | The Sorghum bicolor genome and the diversification of grasses[END_REF]) (see Comparison of assembly protocols in Materials and Methods). We chose these two species because they presented contrasting characteristics (e.g.: eudicot vs monocot,

Accepted Article

This article is protected by copyright. All rights reserved.

highly heterozygote outcrosser vs mainly homozygote selfer, GC-poor vs GC-rich), which allowed testing whether our protocol could be efficient for a wide genome diversity. Especially, highly heterozygote genomes, such as grapevine, should be more challenging than highly homozygote ones, and it is important to know whether protocols should be adapted or not to the expected heterozygosity level. [START_REF] Cahais | Reference-free transcriptome assembly in non-model animals from next-generation sequencing data[END_REF] used reads from multiple individuals to increase the number of assembled contigs as some genes can be lowly or not expressed in some individuals. However, to evaluate the trade-off between maximizing the number of transcripts that can be assembled and optimizing the quality of the assembly by minimizing the complexity of the pool of reads to be assembled, we also compared assembly protocols using reads from multiple accessions or genotypes from the same species (All Individuals) versus protocols using reads from only one individual (One Individual) (Table 2). For the One Individual assembly, we chose the individual with the largest number of reads e.g. EC4

for Sorghum bicolor and VC2 for Vitis vinifera.

Quantitative assessment of the assemblies

For Sorghum bicolor, the number of transcripts predicted from the genome sequence (29,448), was equivalent to the number of contigs produced by the assembly with Abyss One Individual EC4 (27,733), and equivalent to Abyss One Individual EC4 combined with one or two steps of assembly with Cap3 (26,686 and 26,591 respectively). Using reads from all the individuals and an assembly with Abyss (A), Abyss/Cap3 (A/C), or Abyss/Cap3/Cap3 (A/C/C), the number of contigs roughly doubled (62,341, 58,441 and 54,974 respectively). With Trinity EC4, the number of contigs tended towards 70,000 contigs (Figure 1A).

For Vitis vinifera, we confirmed that the number of transcripts predicted from the genomic sequence, Genome (29,971), was close to the number of contigs produced by the assembly with Abyss One Individual VC2 (26,685) and Abyss One Individual VC2

Accepted Article

This article is protected by copyright. All rights reserved.

combined with one or two steps of assembly with Cap3 (24,315 and 24,015 respectively). With Trinity VC2, the number of contigs roughly doubled (48,934). Using all the individuals and Abyss, Abyss/Cap3, and Abyss/Cap3/Cap3 the number of contigs increased 3-4 fold (85,420, 87,125, and 80,982 respectively) (Figure 1B).

For both Sorghum bicolor and Vitis vinifera, the N50 and the mean contig length were higher and closer to those predicted from the genomic sequence (Genome) with all the combinations using one individual (EC4 and VC2) compared to the combinations using all the individuals (Figures 1C and1D). No significant difference in N50 and mean contig length was detected between assembly with Trinity and the different combinations with Abyss using one individual (Figures 1C and1D).

The percentage of contigs mapped on the genome was around 80-90% for Sorghum bicolor and 90-100% for Vitis vinifera (Figure S1). The percentage of predicted genes from the genome sequences mapped by the contigs was about 60% for the One Individual assemblies, and about 70-80% for the assemblies using Trinity and for the All Individuals assemblies (Figures 2A and2B).

Representativeness of the assemblies

For Sorghum bicolor, 85-90% of the total reads from 10 individuals mapped to the assemblies built with Trinity, A/C EC4, A/C/C EC4, and A/C/C EC_All (Figure 3A). This percentage slightly decreased with A EC4 (77%), A/C EC_All (81%), and substantially decreased with A EC_All (53%).

For Vitis vinifera, 85-90% of the total reads from 12 individuals mapped to the assembly built with Trinity, A/C VC2, and A/C/C VC2 (Figure 3B). This percentage slightly decreased with A VC2 (77%), A/C VC_All (67%), A/C/C VC_All (82%), and substantially decreased with A VC_All (30%).

For both Sorghum bicolor and Vitis vinifera, a higher percentage of mapped reads was obtained when the One Individual assembly (produced using either Trinity, A/C or

Accepted Article

This article is protected by copyright. All rights reserved.

A/C/C) was used as reference compared to mapping carried out against CDS transcript sequences predicted from the genome sequence (Figures 3).

The percentage of reads mapped or correctly mapped to the built transcriptome were quite similar, depending on the different individuals. Indeed, we observed less than an 11% difference between individuals for both Sorghum bicolor and Vitis vinifera (Table S1).

For Sorghum bicolor, around 40% of the contigs from the A, A/C and A/C/C One Individual EC4 assemblies were mapped at an average coverage of 8X by reads from the 10 individuals. This percentage decreased notably for the Trinity (15%), A EC_All (23%), A/C EC_All (26%) and A/C/C EC_All (25%) assemblies. For Vitis vinifera, around 66-67% of the contigs from the A, A/C and A/C/C One Individual VC2 assemblies were mapped at an average coverage of 8X by reads from the 12 individuals. This percentage decreased notably for the Trinity (38%), A VC_All (26%), A/C VC_All (37%) and A/C/C VC_All (39%) assemblies (Table S2).

Categorization of the predicted contigs

We applied the categorization of contigs used in [START_REF] Cahais | Reference-free transcriptome assembly in non-model animals from next-generation sequencing data[END_REF] in order to define the assembly protocol minimizing the probability of having erroneous predictions such as chimeras, unassembled alleles and assembled paralogues, which would mislead gene orthology prediction and SNP detection. This categorization was based on the number and the nature of significant BLAST hits of the contigs against a transcriptome predicted from the reference genome annotation (see Classification of BLAST results in Materials and Methods).

The ratio of non-problematic contigs (Full length+Partial+Fragment) was highest with A/C and A/C/C One Individual (>60% for Sorghum bicolor and Vitis vinifera). A/C/C One Individual gave similar results to A/C One Individual: 62 % of non-problematic contigs for A/C EC4 and A/C/C EC4, 64% for A/C VC2 and A/C/C VC2. We slightly

Accepted Article

This article is protected by copyright. All rights reserved.

improved the percentage of non-problematic contigs by combining Abyss with Cap3 compared to Abyss alone with 60% for A EC4 and A VC2. The percentage of nonproblematic contigs decreased for Trinity and All Individuals: 27% for T EC4 and 32% for T VC2, 38% for A EC_ALL and 30% for A VC_ALL, 41% for A/C EC_ALL and A/C/C EC_ALL, 37% for A/C VC_ALL and A/C/C VC_ALL (Figures 4A and4B). With Trinity's assembly using reads from all individuals, we confirmed that the percentage of non-problematic contigs substantially decreased (7% for Sorghum bicolor and 4 % for Vitis vinifera) (Table S3). To assess the impact of isoforms on our classification of contigs, we filtered out isoforms predicted by the Butterfly module of Trinity. From each group of predicted isoforms, we kept the longest one and then we recalculated the classification of contigs. We found that the percentage of non-problematic contigs is still low for Trinity after filtering isoforms: for Sorghum bicolor All Individuals and One individual (15% and 24% respectively), and for Vitis vinifera All Individuals and One individual (12% and 28 % respectively) (Table S4).

We checked whether filters on RPKM could change our results. We compared the percentage of non-problematic contigs for the T EC4, A/C EC4 and A/C EC_ALL assemblies using different filters on RPKM. The A/C EC4 assembly always gave the best score of non-problematic contigs with values greater than 60% (Figure 5).

Conclusions of the assemblies' comparison and choice of the optimal transcriptome assembly strategy

We found that the N50 and mean length of assemblies were higher and closer to those predicted from the genome sequence with all the One Individual combinations compared to the All Individuals combinations (Figures 1C and1D). The percentage of reads from all the individuals that were mapped or correctly mapped to the One Individual assemblies was greater than the percentage of reads mapped to the All

Accepted Article

This article is protected by copyright. All rights reserved.

Individuals contigs (Figures 3).

It is worth noting that the One Individual assemblies produced transcriptome references that were representative of all the sequenced individuals, probably due to fewer erroneous contigs.

We also found that the One Individual assemblies produced by Abyss combined with one step of Cap3 gave the best ratio of non-problematic contigs (Figures 4). We checked that applying filters on read depth of coverage did not notably affect the evaluation results (Figure 5). Lastly, our conclusion was that the optimal protocol was Abyss combined with one step of Cap3, using reads from only one sequenced individual. This protocol led to 95% (Sorghum bicolor) and 97% (Vitis vinifera) of the contigs being mapped on the corresponding reference genome sequence. We favoured this protocol because our aim was to provide reference transcriptomes that could be useful for further comparative and population genomic analyses, even at the expense of transcriptome exhaustivity. It is also interesting that the same protocol performed the best (according to our criteria) for the two tested genomes despite their contrasting characteristics, especially their different heterozygosity level. This is also another advantage to have a protocol robust to change in genome characteristics.

Nevertheless, for other purposes a more exhaustive and detailed representation of the transcriptome might be preferred. In such cases, the Trinity protocol appears to be more suitable, especially for highly heterozygous species.

Assembly of 26 new transcriptomes

We assembled 26 new transcriptomes of cultivated plants and close relative species (listed in Table 1) using the Abyss-Cap3-One Individual protocol. The N50, mean contig length, and percentages of aligned and annotated contigs (Table 3) were comparable to those of other recently published transcriptomes [START_REF] Bhardwaj | Comprehensive transcriptomic study on horse gram (Macrotyloma uniflorum): De novo assembly, functional characterization and comparative analysis in relation to drought stress[END_REF][START_REF] Kudapa | Comprehensive Transcriptome Assembly of Chickpea (Cicer arietinum L.) Using Sanger and Next Generation Sequencing Platforms: Development and Applications[END_REF][START_REF] Munoz-Merida | De novo assembly and functional annotation of the olive (Olea europaea) transcriptome[END_REF][START_REF] Zhang | Comparative analyses of two Geraniaceae transcriptomes using next-generation sequencing[END_REF]. Using our assembly

Accepted Article

This article is protected by copyright. All rights reserved.

protocol, contig numbers ranged from 19,811 in Herrania nitida (a relative of cocoa) to 59,792 in Phillyrea angustifolia (a relative of olive tree) (average = 32,755). Some of our new assemblies (eggplant, pepper, oil palm, olive tree) already have corresponding reference genomes or transcriptomes published. As an example of comparison, we assessed overlaps (see BBMH comparison of transcriptomes in Materials and Methods) between our eggplant assembly and three types of available data sources: the transcripts predicted from the eggplant draft genome SME_r2.5.1 [START_REF] Hirakawa | Draft genome sequence of eggplant (Solanum melongena L.): the representative solanum species indigenous to the old world[END_REF], a previously published transcriptome assembly from RNA-Seq [START_REF] Yang | Comparative transcriptome analysis of eggplant (Solanum melongena L.) and turkey berry (Solanum torvum Sw.): phylogenomics and disease resistance analysis[END_REF], and with a unigene set produced by cDNA sequencing (Vegmark) [START_REF] Fukuoka | Accumulation, functional annotation, and comparative analysis of expressed sequence tags in eggplant (Solanum melongena L.), the third pole of the genus Solanum species after tomato and potato[END_REF]. We found that around 81% (21,574 for a total of 26,752 non redundant sequences) of our eggplant contigs have Best Blast Mutual Hit (BBMH) to at least one of the queried data sources. A core of 8,820 sequences (33%) were commonly found in the three data sets (Figure 6). The published eggplant RNA-Seq assembly (YanXu & al) showed more BBMH (18,703) compared to the transcripts predicted from the draft genome SME_r2.5.1 (15,707) and to the unigene set (11,865). This is probably correlated to a more complete view of the transcriptome provided by RNA-Seq.

Functional annotations

The contigs of the 26 reference transcriptomes were annotated using BLAST, prot4EST [START_REF] Wasmuth | prot4EST: translating expressed sequence tags from neglected genomes[END_REF], and Blast2GO [START_REF] Gotz | High-throughput functional annotation and data mining with the Blast2GO suite[END_REF]. The 26 transcriptomes had more than 60% of the contigs having BLAST hits, and 24 transcriptomes had more than 70% of aligned contigs. For 21 transcriptomes, more than 50% of the contigs had a BLAST hit against the high quality manually curated annotations of the Swiss-Prot database. In addition, 25 transcriptomes had more than 40% of the contigs annotated with Gene Ontology terms (Table 2).

Accepted Article

This article is protected by copyright. All rights reserved.

As an example of biological insights that may be provided by our resource, we assessed GO term enrichment in the group of nine transcriptomes of the Poaceae family (Triticum monococcum, Eremopyrum bonaepartis, Taeniatherum caputmedusae, Sorghum brachypodum, Sorghum bicolor, Pennisetum glaucum, Pennisetum polystachion, Pennisetum alopecuroides, and Digitaria exilis), using the Fisher's Exact Test function implemented in Blast2GO [START_REF] Gotz | High-throughput functional annotation and data mining with the Blast2GO suite[END_REF]. Each transcriptome of the Poaceae group was compared to all other transcriptomes of the group. We filtered with a FDR (False Discovery Rate) value of 0.05. The transcriptome of Sorghum bicolor showed a large number of over-represented terms (8,761) with a slight overrepresentation of terms associated with metabolism (14%) (Figure S2). The transcriptome of Sorghum brachypodum showed a smaller number of GO terms overrepresented (35) (Table S5) with a clear overrepresentation of terms associated with metabolism and binding (Figure 7). The other Poaceae transcriptomes have no more than 5 over-represented GO terms. Note that in this simple analysis "fragment" contigs count twice (or more) for a GO category. This simple analysis could thus be affected by the different species having different proportions of contigs categories. A more detailed analysis of the corresponding contigs and related annotations might lead to new insights into the evolution and/or expression of some specific gene families in Sorghum brachypodum compared to other Poaceae species.

For the 26 new transcriptomes, all the sequences (FASTA format) and annotations (GFF3 format), together with BLAST, prot4EST, and Blast2GO results are freely available at http://arcad-bioinformatics.southgreen.fr/. These transcriptomes and their annotations represent a significant resource for large scale analyses of relationships between the evolution, diversity and function of gene families.

Accepted Article

This article is protected by copyright. All rights reserved.

Claire Billot, Nabila Yahiaoui, Roberto Bacilieri, Bouchaib Khadari, Michel Boccara, Adéline Barnaud, Jean-Pierre Péros, Jean-Pierre Labouisse and Jean-Louis Pham were involved in one or more processes of samples selection and collection, RNA extraction and Illumina library production, data analysis and manuscript preparation.

Figure legends

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

C D

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Figure 1 .

 1 Figure 1. Number of contigs, N50 and mean contig length obtained by the different

Figure 2 .

 2 Figure 2. Percentage of predicted genes mapped by the contigs obtained by the

Figure 3 .

 3 Figure 3. Percentage of the total reads mapped to the assemblies built with the

Figure 4 .

 4 Figure 4. Ratio of non-problematic contigs obtained by the different assembly

Figure 5 .

 5 Figure 5. Effect of RPKM filters on the ratio of non-problematic predicted contigs. The

Figure 6 .

 6 Figure 6. Number of Best Blast Mutual Hits of our eggplant assembly with transcripts

  showing the list size (number of BBMH) and intersection size repartition are located underneath the diagram.

Figure 7 .

 7 Figure 7. Main GO categories of the over-represented GO terms for Sorghum

Figure 1 .

 1 Figure 1. Number of contigs, N50 and mean contig length obtained by the different assembly protocols. The number of contigs increases with Trinity alone and All Individuals (A. Sorghum bicolor, B. Vitis vinifera), the N50 and mean length decreased with All Individuals (C. Sorghum bicolor, D. Vitis vinifera). T, Trinity; A, Abyss; C, Cap3; Genome, transcripts predicted from the genome sequence; EC4, one individual Sorghum bicolor EC4; VC2,one individual Vitis vinifera VC2; EC_ALL, 10 individuals Sorghum bicolor; VC_ALL, 12 individuals Vitis vinifera.

Figure 2 .

 2 Figure 2. Percentage of predicted genes mapped by the contigs obtained by the different assembly protocols. The percentage of predicted genes from the Sorghum bicolor and Vitis vinifera genome sequences that showed an overlap with the contigs was about 60% for One Individual, and 70-80% for Trinity alone and All Individuals: (A) Sorghum bicolor and (B) Vitis vinifera. T, Trinity; A, Abyss; C, Cap3; EC4, one individual Sorghum bicolor EC4; VC2,one individual Vitis vinifera VC2; EC_ALL, 10 individuals Sorghum bicolor; VC_ALL, 12 individuals Vitis vinifera.

Figure 3 .

 3 Figure 3. Percentage of the total reads mapped to the assemblies built with the different protocols. The assemblies with One Individual were the most representative of all the individuals: (A) Sorghum bicolor and (B) Vitis vinifera. T, Trinity; A, Abyss; C, Cap3; Genome, transcripts predicted from the genome sequence; EC4, one individual Sorghum bicolor EC4; VC2,one individual Vitis vinifera VC2; EC_ALL, 10 individuals Sorghum bicolor; VC_ALL, 12 individuals Vitis vinifera.

Figure 4 .

 4 Figure 4. Ratio of non-problematic contigs obtained by the different assembly protocols. The ratio of non-problematic predicted contigs (Full length+Partial+Fragment) was better with AC and ACC One Individual: (A) Sorghum bicolor and (B) Vitis vinifera. T, Trinity; A, Abyss; C, Cap3; EC4, one individual Sorghum bicolor EC4; VC2,one individual Vitis vinifera VC2; EC_ALL, 10 individuals Sorghum bicolor; VC_ALL, 12 individuals Vitis vinifera.

Figure 5 .

 5 Figure 5. Effect of RPKM filters on the ratio of non-problematic predicted contigs. The transcriptome assembly using Abyss with one step of Cap3 (A/C) and reads from the EC4 Sorghum bicolor individual (A/C EC4) maintained the best ratio of non-problematic predicted contigs (% of contigs from the categories Full length+Partial+Fragment) using different filters on RPKM : 0,1,2,3,4 (Reads per Kilobase per Million Reads). T, Trinity; A, Abyss; C, Cap3; EC4, one individual Sorghum bicolor EC4; EC_ALL, 10 individuals Sorghum bicolor.

  1 (Kazusa) was downloaded from ftp://ftp.kazusa.or.jp/pub/eggplant/SME_r2.5.1_cds.fa.gz, the published RNA-Seq transcriptome assembly (YanXu & al) from http://www.ncbi.nlm.nih.gov/Traces/wgs/?download=GBEF01.1.fsa_nt.gz, and the unigene set (Vegmark) from http://cse.naro.affrc.go.jp/fukuoka/EggEST/EggUnigene_16245.fasta.zip.

Table 1 . Summary of sequenced individuals of cultivated plants and close relative species

 1

	Family	Taxon	Crop nam	Life span	Matin g	Number of	Close relati	Accessi on	Germplasm source	Sources of RNA 3
			e		syste	chromo	ve
					m	somes	crop
	Arecace	Elaeis	Oil	Pere	Alloga	2n=32		LM1829	INRAB,
	ae	guineen sis	palm	nnial	my			3-Eg02a	Centre de Recherches Agricoles Plantes Pérennes, Station	80% infloresc ence, 20% leaves
									Pobè, Bénin
		Mauriti		Pere	Alloga	2n=30	Oil	Wild	Field
		a flexuosa		nnial	my		palm	individu al (male), nr. Kourou, French	collection, French Guiana	80% infloresc ence, 20% leaves
								Guiana
	Dioscore	Dioscore	Yam Annu	Alloga	2n=40		CR4885 IRD	80%
	aceae	a		al	my				Montpellier,	infloresc
		rotunda							France	ence,
		ta								20%
										leaves
		Dioscore		Annu	Alloga	2n=40	Yam		CIRAD	80%
		a alata		al	my				collection,	infloresc
									Guadeloupe	ence,
									, French	20%
									West Indies	leaves
		Dioscore		Annu	Alloga	2n=40	Yam		CRB-PT
		a trifida		al	my				(Centre des
									Ressources	80%
									Biologiques_	infloresc
									Plantes	ence,
									Tropicales),	20%
									Guadeloupe	leaves
									, French
									West Indies
	Fabacea	Medica	Alfalf	Pere	Alloga	2n = 4x		Gongno	INRA
	e	go	a	nnial	my	= 32		nd	collection	80%
		sativa							Biological	infloresc
									Ressources	ence,
									Center	20%
									Medicago	leaves
									Montpellier
		Medica		Pere	Alloga	2n = 16	Alfalf	"carnon	INRA	80%
		go		nnial	my		a	"	collection	infloresc

This article is protected by copyright. All rights reserved. This article is protected by copyright. All rights reserved.

There is no clear information regarding the reproductive behaviour of this species. It is quite likely that it is autogamous but with a significant proportion of allogamy (even some cultivated and wild sorghum display significant proportions of allogamy, up to

20% in some cases for cultivated sorghum) 2 Cultivated sorghum is 2n=20 3 Leave means young leave, fruit means immature fruit, and grain means immature grain4 The grape biological materials were provided by the Center for Grape Genetic Resources, at the INRA "Domain de Vassal", Marseillan Plage, France: http://bioweb.ensam.inra.fr/collections_vigne/Home.php?l=EN

The sorghum biological materials were provided by the the CIRAD, Laboratoire de Semences et de Ressources Génétiques, France: LSRG, http://golo.cirad.fr 2 The grape biological materials were provided by the Center for Grape Genetic Resources, at the INRA "Domain de Vassal", Marseillan Plage, France: http://bioweb.ensam.inra.fr/collections_vigne/Home.php?l=EN

Acknowledgements

This study was conducted in the framework of ARCAD (http://www.arcad-project.org), a project funded by Agropolis Fondation under the reference ID ARCAD 0900-001. GS, FH, and SP salaries were supported by ARCAD budget. SC was supported by INRA (Division of Plant Biology and Breeding). The main partner for the sequencing was the Montpellier Genomix regional platform monitored by the CNRS (http://www.mgx.cnrs.fr). We thank Gaëtan Droc for help in setting up the Arcad Bioinformatics Portal. Computational resources were provided by the South Green bioinformatics platform. We thank Dr. Guillaume Besnard and Dr. Pierre Saumitou-Laprade for providing us genetic material of the olive tree out-groups. For yam material and logistics we gratefully acknowledge AfricaRice, Cotonou, and The University of Parakou, Parakou (Benin). We also thank Sally Norton from the Australian Tropical Crops and Forages Collection (QDPI) who provided us with seeds of S. brachypodum. For palm material, RNA extraction and logistics we gratefully acknowledge INRAB Pobè (Benin), PalmElit SAS, H. Adam, J.C. Pintaud, S. Jouannic (IRD), F. Richaud, D. Cros (Cirad) and I. Scotti (INRA). We thank CRC (Cocoa Research Centre) for

providing us with some genetic material (leaves and buds) from T. cacao relatives. We thank Vanuatu Agricultural Research and Technical Centre (VARTC) and Centre de Ressources Biologiques (CRB) for providing us with Dioscorea outgroups. We also thank the institutions that maintain plant/seed collections from which our plant materials originate (listed in the Table1). The technical help of Céline Cardi and Olivier Fouet is acknowledged. We are grateful to Peter Biggins for help in English editing.

Data Accessibility

Raw data were deposited in the NCBI Sequence Read Archive (SRA) under BioProject PRJNA326055.

Sequences assemblies (FASTA), annotations (GFF, Blast XML, Blast2GO .annot, FASTA prot4EST translations): http://arcad-bioinformatics.southgreen.fr.

Scripts used for the analyses: https://github.com/SouthGreenPlatform/arcad-hts.

Author Contributions

Sylvain Glémin, Jacques David, Gautier Sarah and Manuel Ruiz designed research.

Gautier Sarah, Felix Homa, Stéphanie Pointet, Sandy Contreras, François Sabot and Benoit Nabholz developed scripts and performed bioinformatic analyses. Sylvain Glémin, Gautier Sarah and Manuel Ruiz wrote the paper. Sylvain Santoni, François Sabot, Laure Sauné, Morgane Ardisson, Nathalie Chantret, Christopher Sauvage, James Tregear, Cyril Jourda, David Pot, Yves Vigouroux, Hana Chair, Nora Scarcelli,

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.