Caroline Lebaron

Aurélie Rosado

Christopher Sauvage

Camille Gauffier

Sylvie S German-Retana

Benoît Moury

J.-L Gallois

Jean-Luc Gallois
email: jean-luc.gallois@inra.fr

Jean-Luc Gallois

Jean-Luc Gallois

Jean-Luc Gallois

A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability

Keywords: allele mining, RNAseq, lycopersicum spp., eIF4E, Potato virus Y, Tobacco etch virus

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

potential resistance alleles to potyviruses. By analyzing polymorphism among RNAseq data obtained for 20 tomato accessions, ten belonging to the cultivated type Solanum lycopersicum and 10 belonging to the closest related wild species Solanum pimpinellifolium, we isolated one new eIF4E1 allele, in the Solanum pimpinellifolium LA0411 accession, which encodes a potential new resistance allele, mainly due to a polymorphism associated with an amino acid change within eIF4E1 region II.

We confirmed that this new allele, pot1 2 , is indeed associated with resistance to PVY, although with a restricted resistance spectrum and a very low durability potential. This suggests that mutations occurring in eIF4E region II only may not be sufficient to provide efficient and durable resistance in plants. However, our study emphasizes the opportunity brought by RNAseq data to mine for new resistance alleles. Moreover, this approach could be extended to seek for putative new resistance alleles by screening for variant forms of susceptibility genes encoding plant host proteins known to interact with viral proteins.

Introduction

One of the main challenges of plant breeders is to identify new sources of resistance to pathogens.

Crop wild relatives are the main source of such resistances that can be monitored through the direct phenotyping of collection of accessions with the pathogen of interest. However this approach is both time-consuming and costly. Allele mining, which entails looking for polymorphism within genes of interest, represents a complementary approach to plant phenotyping in order to isolate new sources of resistance to pathogens among the natural variation. Viruses are characterized by a particularly small genome and rely on host factors to infect their host (Fraser, 1990). Hence, polymorphism within those host factors can be associated with the plant resistance, a mechanism dubbed as lossof-susceptibility (Pavan et al., 2010;van Schie & Takken, 2014).

Such approach has been exemplified by the characterization of resistance to viruses based on the eukaryotic translation initiation factors eIF4E (Robaglia & Caranta, 2006;Wang & Krishnaswamy, 2012). Translation initiation factors 4E (i.e. eIF4E and the isoform eIFiso4E) are essential components of the eukaryotic cell that initiate translation by binding to the mRNA cap structure at the 5' end of most mRNAs. They are encoded by a small multigene family [START_REF] Browning | Mechanism of cytoplasmic mRNA translation[END_REF]. In the last decade, those factors have been shown to be associated with resistance to a broad range of single-stranded positive-sense RNA viruses and especially potyviruses, including Potato virus Y (PVY) and Tobacco etch virus (TEV) (Wang and Krishnaswamy, 2012). Natural resistance to potyviruses in lettuce [START_REF] Nicaise | The eukaryotic translation initiation factor 4E controls lettuce susceptibility to the Potyvirus Lettuce mosaic virus[END_REF], pepper (Ruffel et al., 2002), pea (Gao et al., 2004) and tomato (Ruffel et al., 2005) has been shown to rely on non-synonymous substitutions in the eIF4E coding sequences.

Those substitutions are mainly located in two regions, named I and II, located near the cap binding pocket in the eIF4E (Robaglia & Caranta, 2006). However, the precise involvement of those two regions in resistance remains unclear. In pepper, a large set of eIF4E1 resistance alleles has been characterized (see below). While mutations in region I of pepper eIF4E1 are associated with resistance to PVY, the additional mutations in region II might be associated with enlargement of the Comment citer ce document : Lebaron, C., Rosado, A., Sauvage, C., Gauffier, C., German-Retana, S., Moury et al., 2008;Yeam et al., 2007). However, no natural allele harboring only resistance-associated mutations in region II has been isolated so far among the pepper natural variation. Beside the need for new resistance alleles, it is also crucial to determine the minimal set of mutations within eIF4E allowing the broadest resistance spectrum possible, in the wake of new breeding technologies such as CRISPR/Cas9 [START_REF] Andersen | Feasibility of new breeding techniques for organic farming[END_REF]. So far, the proof-ofconcept of using this technology to develop resistance has been illustrated by knocking out eIF4E and eIFiso4E genes in Cucumis sativus and Arabidopsis thaliana, respectively (Chandrasekaran et al., 2016;Pyott et al., 2016). However, recent work in tomato suggests that, given the redundancy effect among eIF4E genes, CRISPR/Cas9 technology should rather be used to design functional alleles by introducing non-synonymous mutations in the eIF4E coding sequence, rather than using null alleles [START_REF] Lebaron | A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability[END_REF].

Given their central role in resistance to RNA-viruses, eIF4E family members have also been a target of choice for proof-of-concept studies as well as for a large number of allele mining strategies. In pepper, direct sequencing of eIF4E1 cDNAs sampled on different accessions has characterized 9 pvr2 recessive resistance alleles to PVY and TEV (Charron et al., 2008;Ruffel et al., 2002). Similar approaches using the EcoTILLING strategy or the High Resolution Melting (HRM) analysis have extended this repertoire up to 22 alleles (Ibiza et al., 2010;Jeong et al., 2011). Such studies have also been carried out in other crops such as melon [START_REF] Nieto | EcoTILLING for the identification of allelic variants of melon eIF4E, a factor that controls virus susceptibility[END_REF], barley (Hofinger et al., 2011) or pea [START_REF] Konecna | Geographical gradient of the eIF4E alleles conferring resistance to potyviruses in pea (Pisum) germplasm[END_REF], demonstrating the usefulness of the approach. In the last decade, the development of high-throughput sequencing, especially RNAseq, has provided large amounts of publicly available data and precise knowledge of the gene expression patterns in major crop plants, as well as giving access to the polymorphism associated with amino acid (AA) changes in the encoded proteins. Those data represent a largely untapped resource to mine for new alleles of interest [START_REF] Barabaschi | Next generation breeding[END_REF]. However, it is important to check whether the available data often collected to provide insight on plant evolution and analysis of speciation, in tomato for example (Koenig et al., 2013;Pease et al., 2016), are well suited for the specific purpose of allele mining. (hereafter PI24). This allele is associated with a large resistance spectrum to several isolates of PVY and TEV [START_REF] Legnani | Evaluation and inheritance of the Lycopersicon hirsutum resistance against potato virus Y[END_REF][START_REF] Legnani | Identification and characterization of resistance to tobacco etch virus in Lycopersicon species[END_REF]. Because PI24 is an accession of the wild species Solanum habrochaites, introgression of the resistance into cultivated tomato is difficult [START_REF] Bernacchi | An interspecific backcross of Lycopersicon esculentum × L. hirsutum: linkage analysis and a QTL study of sexual compatibility factors and floral traits[END_REF]Gauffier et al., manuscript in preparation). In this study, we investigated

Results

Analysis of RNAseq data among 20 tomato accessions allows the characterization of a new eIF4E1 allele in the S pimpinellifolium LA0411 accession.

In plants, translation initiation factors are encoded by a small multigene family. In tomato, we previously identified three genes encoding initiation factors 4E, including one isoform eIFiso4E and two eIF4E homologs eIF4E1 and eIF4E2 that are involved in susceptibility to potyviruses (Piron et al., 2010;Ruffel et al., 2005). Verification of the eIF4E2 annotation on the Solgenomics database (https://solgenomics.net/) revealed that the exon structure of the eIF4E2 gene differs from the gene characterized by the substitution of two of the eight tryptophan residues conserved in eIF4E proteins by phenylalanine and tyrosine, respectively (Ruud et al., 1998;Supplemental Figure S1).

Recently, transcriptome data of a set of 20 tomato accessions was obtained through the Arcad project (http://www.arcad-project.org/) using paired-end HiSeq2000 sequencing, aiming at deciphering the molecular footprints of domestication in crop species, such as tomato or rice. The experimental design was based on a comparative approach of 10 accessions belonging to cultivated tomato Solanum lycopersicum and S.lycopersicum cerasiforme and 10 accessions from the closest wild species Solanum pimpinellifolium (Supplemental Table 1). The cDNA sequences for the 4 genes encoding initiation factors, namely eIF4E1, eIF4E2, eIFiso4E and nCBP, were retrieved (Table 1) and polymorphisms called according to the approach detailed in [START_REF] Nabholz | Transcriptome population genomics reveals severe bottleneck and domestication cost in the African rice (Oryza glaberrima)[END_REF] were sought after.

Nucleotide sequences could be retrieved for the four genes in the 20 accessions, showing that those genes are well expressed in leaf tissues used to create the cDNA libraries. No non-synonymous amino-acid changing substitution was found among eIF4E2, eIFiso4E and nCBP coding sequences.

Moreover, no nucleotide polymorphism at all could be found in either eIF4E2 or eIF4iso4E, whereas two silent substitutions were found in the nCBP coding sequence, including one common to all S. pimpinellifolium species sequenced. In comparison, four new haplotypes within the eIF4E1 coding sequence were characterized for tomato eIF4E1 and, as previously shown for the eIF4E1 homolog in pepper, all polymorphisms were associated with AA changes within the eIF4E1 coding sequence (Figure 1). For haplotypes 1 to 3, only one mutation was found, that was located outside the regions I and II where mutations associated with resistance are usually located (Robaglia & Caranta, 2006).

More significantly, the LA0411 accession (haplotype 4) displayed two non-synonymous substitutions 7 within eIF4E1 sequence, including one mutation located at position 112 in the region II, where an aspartic acid (D) is replaced by a glycine (G) (D112G, Figure 1).

The presence of this polymorphism was confirmed in the LA0411 eIF4E1 allele by Sanger sequencing.

The rest of the study was focused on the LA0411 allele, as it encodes a protein with an AA change within the region II but not in region I. It was of particular interest to see whether the region II mutation by itself provides resistance to potyviruses.

The eIF4E1 allele of accession LA0411 is a new pot1 2 resistance allele to two PVY isolates

We first checked whether the S. pimpinellifolium LA0411 accession displayed any resistance to potyviruses. Thus, LA0411 plants were mechanically inoculated with three isolates of PVY (LYE90, N605 and SON41g) and three isolates of TEV (CAA10, S103 and HAT). The resistance phenotype was assayed by titrating viral accumulation by ELISA at 21 days post-inoculation, in the systemicallyinfected leaves. Susceptible S. lycopersicum plants M82, as well as the highly resistant S.

habrochaites PI24 accession, which harbors the eIF4E1-pot1 allele, were used as positive and negative controls, respectively. Our results showed that the LA0411 plants display resistance to the PVY isolates LYE90 and N605 and is susceptible to PVY SON41g as well as to all three TEV isolates tested (Table 2 and supplemental figure 2). Therefore, LA0411 is associated with a reduced but significant resistance spectrum to potyviruses.

The genetic basis of LA0411 resistance to potyviruses was further assessed using PVY N605 isolate.

LA0411 plants were manually crossed with the susceptible S. lycopersicum M82 accession. F1 plants were crossed to each parent to produce resistant and susceptible backcrosses, respectively. F2 were produced by selfing F1 plants. All populations were tested for resistance or susceptibility to PVY N605 (Table 3). Both the F1 and the susceptible backcross progeny plants were found to be fully susceptible to PVY N605, suggesting the presence of a recessive resistance gene. In comparison, Finally, an allelism test was carried out between the broad resistance allele eIF4E1-pot1 that originates from the wild accession PI24 (Parrella et al., 2002) and the LA0411 eIF4E1 allele. Wild species belonging to S. habrochaites, such as PI24, are notably difficult to cross with cultivated tomato [START_REF] Bernacchi | An interspecific backcross of Lycopersicon esculentum × L. hirsutum: linkage analysis and a QTL study of sexual compatibility factors and floral traits[END_REF]. Therefore, near-isogenic lines that have been generated by introgressing the eIF4E1-pot1 allele into the elite Mospomorist (S. lycopersicum) cultivar were used, hereafter named NIL-pot1 (Gauffier et al., in preparation). F1 plants issued from the cross between NIL-pot1 and LA0411 were all resistant to PVY N605 (n=20 plants) suggesting that the monogenic recessive resistance in LA0411 is allelic to eIF4E1-pot1, and that the resistance in the S.

pimpinellifolium LA0411 is caused by mutations in its eIF4E1 allele. Consequently, we propose to name this allele as eIF4E1-pot1 2 . To confirm this allelism test, F1 plants generated between NIL-pot1 and LA0411 were self-crossed and F2 plants were tested for their resistance to PVY-N605. Although the F2 plants were mostly resistant to PVY-N605 as expected given the allelism between eIF4E1-pot1 and eIF4E1-pot1 2 (38 resistant plants out of 40 plants inoculated with PVY-N605), two plants (representing 5%) were found to be fully susceptible to the virus. This low occurrence of susceptibility could be caused by resistance breakdown of the eIF4E1 resistance allele. Consequently, we further investigated the durability of eIF4E1-pot1 2 mediated resistance to PVY-N605.

The resistance associated with eIF4E1-pot1 2 is easily overcome by PVY N605-derived variants

Partial breakdown of the resistance harbored by the S. pimpinellifolium accession LA0411 could have been missed due to the small numbers of plants assayed. Therefore, two more resistance tests were carried out, by assaying at least 50 plants per genotype. All M82 plants were found to be susceptible while all PI24 plants were resistant to the virus. In comparison, LA0411 plants showed a surprising high level of susceptibility ranging from 43 % to 56% of the plants (Table 4). et al., 2008;[START_REF] Moury | Mutations in Potato virus Y genome-linked protein determine virulence toward recessive resistances in Capsicum annuum and Lycopersicon hirsutum[END_REF]. Sequencing of the RT-PCR-amplified VPg-cistron of the PVY N605 isolate after its propagation in three independent susceptible M82 plants did not reveal any mutation in comparison with the original inoculum. Then, 25 independent LA0411 plants infected by PVY-N605

were sampled and the corresponding VPg cistrons sequenced to identify potential mutations associated with the gain of virulence. All sequences showed polymorphisms associated with AA changes within the VPg coding region. These mutations were a substitution of a leucine by a serine at position 115 in all of the sequenced progenies (L115S) and an additional substitution of an isoleucine by a valine at position 139 (I139V) in 9 of the sequenced progenies. This result confirms that resistance-breaking variants derived from PVY N605 were detected after a first passage on LA0411 plants. Notably, the mutation affecting L115 had previously been characterized as a mutation consistently associated with resistance-breaking in pepper for the related, PVY SON41p isolate [START_REF] Ayme | Different mutations in the genome-linked protein VPg of potato virus Y confer virulence on the pvr2(3) resistance in pepper[END_REF].

To further confirm that LA0411 susceptibility is caused by resistance-breaking PVY N605 variants that emerged in the progeny, back-inoculations were carried out. Three kinds of viral inoculum were used: a PVY N605 viral progeny propagated on susceptible M82 plants, compared with two PVY N605 variant progenies propagated in LA0411 plants and carrying the L115S and L115S+I139V substitutions, respectively. All inoculated M82 plants were susceptible to the three isolates, whereas all PI24 inoculated plants were resistant (Table 5). The later result shows that the mutations acquired during the first passage on LA0411 did not confer to the PVY N605 variants the ability to overcome eIF4E1-pot1 mediated-resistance in the broadly resistant line PI24. When the evolved viral isolates present in the LA0411 symptomatic plants after a first passage were back-inoculated onto LA0411 plants, 93% of the plants were infected. Therefore, the second passage allowed the increase in the number of infected plants (ranging from 62% to 93%) confirming that the L115S and L115S+I139V

variants were breaking the resistance associated with eIF4E1-pot1 2 . The remaining non-infected Altogether, these results show that the eIF4E1-pot1 2 allele is associated with a very low durability to PVY-N605.

The PVY-N605 resistance-breaking isolates use both eIF4E1 and eIF4E2 in tomato

We showed previously that although the natural allele eIF4E1-pot1 was associated with resistance to most PVY and TEV isolates, a null TILLING allele knocking out (KO) eIF4E1 was associated with a very narrow resistance spectrum (Piron et al., 2010). We further showed that resistance could be restored, including to PVY-N605, by combining null mutations affecting both eIF4E1 and eIF4E2, 11 the eIF4E1-pot1² allele with a null eIF4E2 TILLING allele will be carried out to confirm this hypothesis:

the resulting genotype is expected to be susceptible to the RB PVY-N605 isolates.

Discussion

In the present study, by using RNAseq data collected from twenty tomato accessions, including 10 from the wild related species Solanum pimpinellifolium, we could assess for the first time the sequence variability among all eukaryotic translation initiation factors 4E, including the isoform eIFiso4E and the atypical related protein nCBP that had not been described yet in tomato. This new set of data is crucial in the light of the comparison with data already put forward in pepper (Capsicum spp.) for eIF4E1 (Cavatorta et al., 2008;Charron et al., 2008;Ibiza et al., 2010). Indeed, both pepper and tomato are close relatives among the Solanaceae family and both are host of PVY and TEV potyviruses. In both Capsicum and Solanum species, recessive resistance to these potyviruses has been characterized as relying on natural variants of eIF4E1. Finally, both plant species encode the same set of eukaryotic translation initiation factors 4E: eIF4E1, a closely related eIF4E2 gene, the eIFiso4E isoform and the nCBP [START_REF] Lebaron | A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability[END_REF] this study; Gallois and Caranta, unpublished results).

Our results suggest contrasted evolution patterns within the eukaryotic translation initiation factor 4E gene family. As shown previously in pepper, the eIF4E1 gene displayed the largest number of nucleotide polymorphisms associated with AA changes. Overall, the screening of the 20 tomato accessions from the ARCAD data set allowed us to find polymorphism within 7 accessions and identify 4 new haplotypes coding for 4 distinct new eIF4E1 proteins. Those polymorphisms affected five positions within the eIF4E1 coding sequence and were all associated with AA changes. These results are in accordance with the analysis of polymorphism of the Capsicum eIF4E1 (Cavatorta et al., 2008;Charron et al., 2008;Ibiza et al., 2010) et al., 2009;Quenouille et al., 2016;Ruffel et al., 2006). Therefore, this nucleotide variability in the pepper eIFiso4E could be driven by an involvement in susceptibility to potyviruses, and on the contrary, eIFiso4E would not be involved in susceptibility to potyviruses in tomato. Indeed all eIFiso4E loss-of-function tomato mutants studied so far remain susceptible to potyvirus infection (Piron et al., 2010). Finally, we report for the first time an analysis of nucleotide diversity within the eIF4E2 coding sequence. No polymorphism could be detected among the eIF4E2 coding sequence for any of the 20 tomato accessions tested. eIF4E2 has been shown to be a susceptibility factor for potyviruses, but only when the eIF4E1 gene was knocked out, through RNAi or EMS-induced mutation [START_REF] Lebaron | A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability[END_REF][START_REF] Mazier | Knock-down of both eIF4E1 and eIF4E2 genes confers broad-spectrum resistance against potyviruses in tomato[END_REF]. More precisely, comparing natural functional eIF4E1 variants and non-functional KO TILLING mutants in tomato allowed us to show that the presence of a eIF4E1 natural resistance allele makes the eIF4E2 factor unavailable for the potyvirus. Consequently, the stability of eIF4E2 coding sequence is fully consistent with an absence of selection pressure from the viruses on eIF4E2.

Similar results were obtained on the pepper eIF4E2 sequence variability (Charron, 2007), suggesting that a similar mechanism of redundancy between eIF4E1 and eIF4E2 could be found in this related species.

This study allowed us to characterize in Solanum pimpinellifolium accession LA0411 one new eIF4E1 allele whose sequence suggests that it could provide a new resistance allele, due to the presence of 2006). This mutation is also present on its own in two cultivated accessions from our set (LA1420 and LA0409). Both LA1420 and LA0409 plants were found to be fully susceptible to PVY N605 (10 susceptible plants out of 10) suggesting that this AA change does not play a role in resistance (data not shown). The second mutation, D112G, is located in region II, a short region predicted to locate near the cap-binding pocket in the eIF4E 3D structure and involves a substitution of an acidic residue by a non-polar one. Because it originates from a close wild relative to cultivated tomato, this allele can be easily introgressed into the cultivated tomato. However, compared with the eIF4E1-pot1 resistance allele, which originates from the PI24 Solanum habrochaites accession and confers resistance to PVY and TEV, eIF4E1-pot1 2 is associated with a narrow resistance spectrum as it confers PVY resistance with limited durability and no resistance to TEV. These results suggest also that efficient resistance to potyviruses in Solanaceae may require simultaneous AA changes in two distant regions of the eIF4E1 protein (regions I and II): these results must be taken into account for strategies involving new breeding techniques such as CRISPR/Cas9 to engineer resistance alleles that may require two rounds of CRISPR/Cas9 gene modifications. Alternatively, one may consider combining the use of natural variation and genome modification by using the eIF4E1-pot1 2 allele as a template to introduce mutations in region I by CRISPR/Cas9.

The development and evolution of high-throughput sequencing technologies, such as RNAseq, have a tremendous impact on plant breeding by allowing access to the genome-wide genetic diversity, molecular marker development and genome-wide association studies at a reasonable cost [START_REF] Barabaschi | Next generation breeding[END_REF]Goodwin et al., 2016). Our study demonstrates also the straightforward use of RNAseq data to mine for alleles of interest, and can be implemented using new sets of data, on a et al., 2016). However, the limitations of next generation sequencing data should be kept in mind, including the presence of false negative single nucleotide polymorphisms (i.e. the lack of detection of actual polymorphisms), a problem already acknowledged for clinical genetic diagnostic (Huang et al., 2015;Park et al., 2015). Such false negative rates have been assessed as high as 4% by comparing RNAseq data with Sanger sequences in Brassica rapa and could account for missing potentially interesting polymorphisms (Devisetty et al., 2014). Therefore, the lack of detected polymorphism does not necessarily reflect a wild-type allele. For example, Sanger sequencing confirmed that the Cervil cultivated tomato accession, included in our panel, contains an eIF4E1 variant (D112N) that was not detected in our RNAseq data screening.

Overall, the use of RNAseq data potentially allows the simultaneous analysis of all members from the same gene family. That is of particular interest for plant translation initiation factors 4E, given the large redundancy among the genes encoding those factors in plant genome, as well as the specific recruitment of the different eIF4E isoforms by various species of potyviruses as exemplified in Arabidopsis thaliana (Duprat et al., 2002;[START_REF] Nicaise | Coordinated and selective recruitment of eIF4E and eIF4G factors for potyvirus infection in Arabidopsis thaliana[END_REF]Sato et al., 2005). Furthermore, such studies could be extended by investigating a large number of susceptibility candidate genes, such as the chloroplastic phosphoglycerate kinase 2 (cPGK2) that has been shown to be associated with resistance to the potyvirus Watermelon mosaic virus in Arabidopsis thaliana, or the protein disulfide isomerase Like 5-1 (HvPDIL5-1), a gene involved with resistance to bymoviruses in barley [START_REF] Ouibrahim | Cloning of the Arabidopsis rwm1 gene for resistance to Watermelon mosaic virus points to a new function for natural virus resistance genes[END_REF]Yang et al., 2014). Because of their limited genome, viruses recruit host factors to complete all the steps of their infectious cycle and are therefore engaged in many interactions with host proteins. Several studies have allowed the characterization of plant proteins interacting with viral proteins, mainly using yeast two hybrid protein-protein interaction systems (for review see

Elena & Rodrigo 2012) or affinity purification [START_REF] Martinez | Interaction network of tobacco etch potyvirus NIa protein with the host proteome during infection[END_REF]. Those host genes are potentially encoding susceptibility factors, and variants of those genes could possibly encode new resistance factors, mirroring the example of eIF4E proteins in potyvirus resistance. This reverse genetics 15 approach could allow extending the portfolio of resistance genes based on loss-of-susceptibility and could define potentially resistant accessions that could then be used as sources for transfer to crops.

Methods.

Plant materials M82 (S. lycopersicum) was used as the susceptible control. The resistant accession S. habrochaites PI247087 (PI24) has been described previously in Parrella et al. (2002). LA0411 accession seeds were provided by INRA CRB (Centre de Ressources Biologiques, http://www6.paca.inra.fr/gafl). The TILLING lines KO for eIF4E1 (∆E1) and eIF4E2 (∆E2) were previously obtained in a S. lycopersicum M82 background (Piron et al., 2010) and the double mutant (∆E1, ∆E2) has been previously described [START_REF] Lebaron | A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability[END_REF].

eIF4E sequence analysis

The RNAseq data were obtained in the framework of the ARCAD project aiming at deciphering the molecular footprints of domestication through a comparative genomic approach. Towards this objective, the experimental design relied on 10 crop and 10 wild accessions from 2 species (S. lycopersicum and S. pimpinellifolium) chosen on the basis of their molecular diversity. A standardized protocol was followed for each of these species in terms of RNA extraction, RNA sequencing and bioinformatic analyses as described in [START_REF] Nabholz | Transcriptome population genomics reveals severe bottleneck and domestication cost in the African rice (Oryza glaberrima)[END_REF]. Among these species, crop tomato accessions and wild accessions were grown in the greenhouse under classical conditions. Briefly RNA was extracted from young leaves, fruits and flowers that were pooled in a 65/20/15 proportion following their quantification on an Agilent Bioanalyser. Once quantified, RNAseq libraries were prepared using TruSeq RNA sample Preparation v2 kit (Illumina Inc., USA) according to the manufacturer's protocol and individually tagged (6bp). These paired-end libraries (2x101bp) were sequenced and the raw reads were checked and cleaned using FastQC (0.10.2, (v1.8.5, Martin, 2011, with the options -q 20 -e 0.1 -O 7 -m 20) tools. Then, these reads were mapped onto the reference sets of ORF corresponding to the v2.4 of the tomato reference genome annotation (n=34,727 ORFs), using BWA (v0.7.3, Li and Durbin, 2009, using Raw reads are hosted online on the Southgreen platform (http://arcad-bioinformatics.southgreen.fr).

For sequence confirmation, Total RNA was isolated from tomato leaves using TRI-reagent (Sigma-Aldrich, Saint Louis, MO, USA). eIF4E1 cDNAs from LA0411 were obtained by RT-PCR with primer pairs F 5'-ATGGCAGCAGCTGAAATGGAGAGA-3'/R 5'-CTATACGGTGTAACGATTCTTGGC-3'. To compare polymorphisms between the accessions, amplification products were sequenced by Genoscreen (Lille, France) and the corresponding protein sequences were aligned using Clustal Omega (http://www.ebi.ac.uk/Tools/msa/clustalo; website) and BoxShade (http://www.ch.embnet.org/software/BOX_form.html; website).

Virus isolates and infection assays

The PVY isolates LYE90v [START_REF] Moury | Mutations in Potato virus Y genome-linked protein determine virulence toward recessive resistances in Capsicum annuum and Lycopersicon hirsutum[END_REF]), N605 (Parrella et al., 2002) and SON41g (Charron et al., 2008) and the TEV isolates CAA10 (Charron et al., 2008), S103 (Ruffel et al., 2005) and HAT (Schaad et al., 2000) were propagated on Nicotiana tabacum cv Xanthi before inoculating 14-day-old tomato plants. The accumulation of PVY and TEV viruses in non-inoculated upper leaves was then assayed 21 days post inoculation by DAS-ELISA using respectively anti-PVY (Sediag, France) and anti-TEV (Sediag, France) antisera and detection sets. Non-inoculated plants were used as healthy controls. Mean and standard errors of absorbance values at 405nm of samples from at least 6 independent plants per

- - S. pimpinellifolium LA0411 0 7 0:1 - - (M82 X LA0411)F1 20 0 1:0 - - BC1S = (M82 X LA0411)F1 X M82 20 0 1:0 - - BC1R = (M82 X LA0411)F1 X LA0411 20 20 1:1 - - (M82 X

  uncovering a redundancy effect between eIF4E1 and eIF4E2[START_REF] Lebaron | A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability[END_REF]. To investigate how the Resistance Breaking (RB) PVY-N605 variants are able to overcome the eIF4E1-pot1 2 mediated resistance, we looked at the resistance status of previously characterized tomato TILLING loss-of-function mutants[START_REF] Lebaron | A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability[END_REF] Piron et al., 2010) towards those RB variants. The reference isolate PVY-N605 and the two RB isolates, respectively harboring L115S and L115+I139V substitutions in the VPg, were propagated on susceptible M82 plants and inoculated on single KO mutant affecting eIF4E1 and eIF4E2 (hereafter 4E1 and 4E2 plants, respectively) and on plants combining both mutations 4E1 4E2 (Table6). 4E1 and 4E2 plants were susceptible to all three PVY isolates tested, showing that PVY-N605 isolate and its derived eIF4E1-pot1 2 resistance breaking PVY-N605 L115S and L115+I139V variants, are able to recruit either eIF4E1 or eIF4E2 to infect tomato. In comparison, plants combining mutations in both eIF4E1 and eIF4E2, 4E1X4E2, in addition to being resistant to PVY-N605, were also fully resistant to the RB variants PVY-N605 L115S and L115+I139V. As suggested for the pepper/ PVY and pepper/ TEV pathosystems, these results suggest that the mutations in the VPg of the two evolved PVY-N605 isolates allow the virus to recruit the mutated form of eIF4E1-pot1² rather than hijacking new host factors in the plant. Combination of Comment citer ce document : Lebaron, C., Rosado, A., Sauvage, C., Gauffier, C., German-Retana, S., Moury, B., Gallois, J.-L. (Auteur de correspondance) (2016). A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability. Journal of General Virology, 97 ( 11), 3063-3072. , DOI : 10.1099/jgv.0.000609 Downloaded from www.microbiologyresearch.org by IP: 147.100.66.219 On: Mon, 03 Oct 2016 12:08:58

  the command line 'aln -e -1 -O 11 -n 2 -m 2000000 -d 10 -E 4 -l 32 -M 3 -k 2 -q 0 -R 30 -t 4 -i 5 -o 1') by allowing no more than two mismatches per read, no multiple hits and the local realignment option. Finally GATK (McKenna et al., 2010) called simultaneously the SNP. The VCFtools (v0.1.10, Danacek et al., 2011) filtered the output variant calling file (vcf) to retain sites showing a minimal coverage of 8x per individual and a MAF of 0.05.

Figure 1 .

 1 Figure legend: Figure 1. New eIF4E1 haplotypes characterized in this study. (a) Sequence alignment of eIF4E1 proteins encoded by S. lycopersicum M82 and S. habrochaites PI24 and PI13 alleles (Ruffel et al., 2005) as well as those encoded by the four new haplotypes characterized in this study. Regions I and II (defined by Robaglia & Caranta 2006) are represented in red and green, respectively. (b)Table of the accession where the four new haplotypes have been characterized.

  ...................................................K...... haplotype2 1 .........................................................V.. haplotype3 1 ............................................................ haplotype4 1 .....................................................K...... PI13 1 ............................................................ PI24 1 ...............................................F................................................................... haplotype2 61 ............................................................ haplotype3 61 ............................................................ haplotype4 61 ...................................................G........ PI13 61 ........S...............L................................... PI24 61 .......KS.......D.......L.......................I...........M82 121 EPKWEDPVCANGGTWKMSFSKGKSDTSWLYTLLAMIGHQFDHGDEICGAVVSVRAKGEKI haplotype1 121 ............................................................ haplotype2 121 ............................................................ haplotype3 121 ..............................M............................. haplotype4 121 ............................................................ PI13 121 ..Q......................................................... PI24 121 ..Q......................................................... M82 181 ALWTKNAANETAQVSIGKQWKQFLDYSDSVGFIFHDDAKRLDRNAKNRYTV haplotype1 181 ................................................... haplotype2 181 ................................................... haplotype3 181 ................................................... haplotype4 181 ................................................... PI13 181 ...........................................S....... PI24 181 ...........................................S.......

Figure 1

 1 Figure 1 Click here to download Figure Figure 1 LEBARON without legenda.docx

  Compared to the large amount of eIF4E resistance alleles discovered in Capsicum Spp., only one resistance allele of the eIF4E1 tomato orthologous gene was discovered in the accession PI247087

Comment citer ce document : Lebaron, C., Rosado, A., Sauvage, C., Gauffier, C., German-Retana, S., Moury, B., Gallois, J.-L. (Auteur de correspondance) (2016). A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability. Journal of General Virology, 97 ( 11), 3063-3072. , DOI : 10.1099/jgv.0.000609 Downloaded from www.microbiologyresearch.org by IP: 147.100.66.219 On: Mon, 03 Oct 2016 12:08:58 5

  Charron et al., 2008), and, consequently, a similar process could occur in tomato. In contrast with the high level of polymorphism among eIF4E1, no polymorphism was detected in eIF4E2 or in eIFiso4E sequences. The results on eIFiso4E nucleotide variability in tomato differ from the relatively high diversity of polymorphism discovered among Capsicum spp.

	accessions, again associated with AA changes (Ibiza et al., 2010). Interestingly, several pepper
	accessions harbor a natural knock-out of the eIFiso4E gene, due to a deletion occurring within the
	gene. This loss-of-function, when combined with a eIF4E1-mediated pvr2 resistance allele is
	associated with resistance to several isolates of Pepper veinal mottle virus and Chilli veinal mottle

and suggest that, as in pepper, positive selection may act on eIF4E1 to fix new alleles. In pepper, this positive selection has been associated with Comment citer ce document : Lebaron, C., Rosado, A., Sauvage, C., Gauffier, C., German-Retana, S., Moury, B., Gallois, J.-L. (Auteur de correspondance) (2016). A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability. Journal of General Virology, 97 ( 11), 3063-3072. , DOI : 10.1099/jgv.0.000609 Downloaded from www.microbiologyresearch.org by IP: 147.100.66.219 On: Mon, 03 Oct 2016 12:08:58 12 coevolution with potyviruses (virus, as well as an increase in pvr2 resistance overcoming (Hwang

  change within region II. Both genetic studies and resistance assays showed that this new allele is indeed a recessive eIF4E1 resistance allele to PVY, hence named pot1 2 . The eIF4E1 protein encoded by this allele is characterized by two AA changes. One of these, V54K, is located at the N-terminal part of the protein, outside regions I and II where mutations associated with resistance to potyviruses are consistently located in eIF4E factors in different plant species(Robaglia & Caranta,

Comment citer ce document : Lebaron, C., Rosado, A., Sauvage, C., Gauffier, C., German-Retana, S., Moury, B., Gallois, J.-L. (Auteur de correspondance) (2016). A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability. Journal of General Virology, 97 ( 11), 3063-3072. , DOI : 10.1099/jgv.0.000609 Downloaded from www.microbiologyresearch.org by IP: 147.100.66.219 On: Mon, 03 Oct 2016 12:08:58 13 an AA

  Comment citer ce document : Lebaron, C., Rosado, A., Sauvage, C., Gauffier, C., German-Retana, S., Moury, B., Gallois, J.-L. (Auteur de correspondance) (2016). A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability. Journal of General Virology, 97 ( 11), 3063-3072. , DOI : 10.1099/jgv.0.000609

	http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and Cutadapt
	16
	Downloaded from www.microbiologyresearch.org by
	IP: 147.100.66.219
	On: Mon, 03 Oct 2016 12:08:58

  Table of the accession where the four new haplotypes have been characterized.

Table 1 . Haplotypes of the tomato translation initiation factors 4E found in the RNAseq ARCAD dataset.

 1

	Tables:
	Gene	Identifier	Chromosome CDS length	Haplotypes	New haplotypes with AA changes
	eIF4E1	Solyc03g005870	3	696	5	4
	eIF4E2	Solyc02g021550	2	663	1	0
	eIFiso4E	Solyc09g090580	9	603	1	0
	nCBP	Solyc10g080660	10	672	3	0
	CDS: Coding DNA Sequence. AA: Amino Acid. The haplotypes references correspond to the sequences from Solanum
	lycopersicum M82 accession .

Table 2 . LA0411 is associated with a narrow resistance spectrum to PVY

 2

	PVY	TEV

Comment citer ce document : Lebaron, C., Rosado, A., Sauvage, C., Gauffier, C., German-Retana, S., Moury, B., Gallois, J.-L. (Auteur de correspondance) (2016). A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability. Journal of General Virology, 97 ( 11), 3063-3072. , DOI : 10.1099/jgv.0.000609 Downloaded from www.microbiologyresearch.org by IP: 147.100.66.219 On: Mon, 03 Oct 2016 12:08:58 22

Table 4 . LA0411 shows lower durability to PVY N605 than PI24

 4

				Accession
		M82		PI24		LA0411
		R	n	R	n	R	n
	Assay 1	0%	57	100%	57	43%	58
	Assay 2	0%	51	100%	60	56%	59
	Results from DAS-ELISA from two independent assays. R represents the percentage of resistant
	plants. n = number of plants assayed. ELISA were performed at 21 dpi.

Table 3 . Inheritance of resistance to PVY N605 in Solanum pimpinellifolium LA0411

 3

	Genotype	Phenotype	Expected ratio	CHI 2	P
		S	R
	S. lycopersicum M82	7	0	1:0

Table 5 . Back inoculation of PVY-N605 evolved isolates on LA0411 plants is associated with a reduction of resistance

 5 ELISA. R represents the percentage of resistant plants. n = number of plants assayed.

	LA0411)F2	26	14	3:1	2.13	0.14

Comment citer ce document : Lebaron, C., Rosado, A., Sauvage, C., Gauffier, C., German-Retana, S., Moury, B., Gallois, J.-L. (Auteur de correspondance) (2016). A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability. Journal of General Virology, 97 ( 11), 3063-3072. , DOI : 10.1099/jgv.0.000609 Downloaded from www.microbiologyresearch.org by IP: 147.100.66.219 On: Mon, 03 Oct 2016 12:08:58 23

Table 6 . Resistance of tomato 4E TLLING mutants to PVY N605 and two RB variants

 6 Resistance (R) to PVY was characterized by DAS-ELISA. For each category, the number of resistant plants out of the total number of plant assayed is given.

			M82 Accession
		4E1	4E2	4E1 4E2
	PVY N605	0 R / 8	0 R / 8	8 R / 8
	L115S	1 R / 7	0 R / 7	7 R / 7
	L115S + I139V	0 R / 8	0 R / 7	8 R / 8

Comment citer ce document : Lebaron, C., Rosado, A., Sauvage, C., Gauffier, C., German-Retana, S., Moury, B., Gallois, J.-L. (Auteur de correspondance) (2016). A new eIF4E1 allele characterized by RNAseq data mining is associated with resistance to PVY in tomato albeit with a low durability. Journal of General Virology, 97 ( 11), 3063-3072. , DOI : 10.1099/jgv.0.000609 Downloaded from www.microbiologyresearch.org by IP: 147.100.66.219 On: Mon, 03 Oct 2016 12:08:58

Acknowledgements

The RNAseq data production was initially supported by Agropolis Foundation under the ARCAD project No 0900-001 (http://www.arcad-project.org/). We thank Renaud Duboscq for additionnal RNAseq data analysis as well as André Moretti and Carole Caranta for providing the NIL-pot1 plants.

We thank the INRA Avignon Centre de Ressources Biologiques (CRB) for providing the LA0411 seeds.

Comment citer ce document : Lebaron, C., Rosado, A., Sauvage, C., Gauffier, C., German-Retana, S., Moury


Alexandre Marchal

Andrés Legarra

Sébastien Tisné

Catherine Carasco-Lacombe

Aurore Manez

Edyana Suryana

INRAB, CRAPP, Pobè Alphonse Omoré

Bruno Nouy

Tristan Durand-Gasselin

Leopoldo Sánchez

Jean-Marc M Bouvet

David Cros
email: david.cros@cirad.fr

P T Socfindo Medan

Benin B Nouy

Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests

Keywords: Elaeis guineensis, genomic selection, multivariate model, empirical data, reciprocal recurrent selection, accuracy

Genomic selection is promising for plant breeding, particularly for perennial crops. Multivariate analysis, which considers several traits jointly, takes advantage of the genetic correlations to increase accuracy.

The aim of this study was to empirically evaluate the potential of a univariate and multivariate genomic mixed model (G-BLUP) compared to the traditional univariate pedigree-based BLUP (T-BLUP) when analyzing progeny tests of oil palm, the world major oil crop.

The dataset comprised 478 crosses between two heterotic groups A and B with 140 and 131 parents, respectively, genotyped with 313 SSR. The traits were bunch number and average bunch weight.

We found that G-BLUP with a genomic matrix based on a similarity index had a higher likelihood than T-BLUP. Also, multivariate G-BLUP improved the accuracy of additive effects (breeding values or general combining abilities, GCAs), in particular for the less heritable trait, and of dominance effects (specific combining abilities, SCAs). The average increase in accuracy was 22.5% for GCAs and 18.7% for SCAs. Using 160 markers in group A and 90 in group B was enough to reach maximum GCA prediction accuracy. The contrasted history of the parental groups likely explained the higher benefit of G-BLUP over T-BLUP for group A than for group B.

Finally, G-BLUP should be used instead of T-BLUP to analyze oil palm progeny tests, with a multivariate approach for correlated traits. G-BLUP will allow breeders to consider SCAs in addition to GCAs when selecting among the progeny-tested parents.

Introduction

Oil palm (Elaeis guineensis Jacq.) is the main oil crop in the world. It bears fruit bunches all year long, and palm oil is extracted from the mesocarp of the fruits. Bunch production is a key component of oil yield, and results from the product of two negatively correlated traits, bunch number (BN) and average bunch weight (ABW) [START_REF] Corley | Selection and breeding[END_REF]. Commercial oil palms are hybrids between two heterotic groups called A and B. Group A is mostly made up of the Deli population (Asia) and group B of various African populations. Group A palms have a few heavy bunches whereas group B palms have many small bunches, resulting in heterosis of bunch yield in A × B hybrids. This led to the choice of a reciprocal recurrent selection (RRS) breeding scheme in the 1950s [START_REF] Gascon | Caractéristique de la production d'Elaeis guineensis (Jacq.) de diverses origines et de leurs croisements. Application à la sélection du palmier à huile[END_REF][START_REF] Meunier | Le schéma général d'amélioration du palmier à huile à l'IRHO[END_REF]. RRS involves progeny tests in which group A and group B parents are crossed to estimate their general combining ability (GCA), i.e. half their breeding value in hybrid crosses, for each yield component, from the phenotype of their hybrid progenies. So far, parental GCAs are obtained using an univariate mixed-model analysis (i.e. considering one trait at a time) taking pedigree information into account [START_REF] Soh | Ranking parents by best linear unbiased prediction (BLUP) breeding values in oil palm[END_REF][START_REF] Purba | Prediction of oil palm (Elaeis guineensis, Jacq.) agronomic performances using the best linear unbiased predictor (BLUP)[END_REF]. The accuracy of the GCAs (i.e. the correlation between the estimated and the true GCAs) is high, reaching around 0.9 for all yield components (Cros et al. 2015b). However, the progeny tests require a long generation interval (around 20 years) and low selection intensity (less than 200 individuals tested per parental group and generation).

Genomic selection (GS) aims to predict genetic values of candidate individuals. In particular, GS can be applied on candidate individuals without data records, by using their genotype with high density molecular markers and a model calibrated with a training set made of individuals with records and marker data [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF]. GS is then particularly promising when traditional breeding requires extensive phenotyping, like progeny tests, as in this case GS makes it possible to reduce the generation interval and to increase selection intensity. The potential of GS is particularly high for perennial crops [START_REF] Grattapaglia | Breeding forest trees by genomic selection: current progress and the way forward[END_REF][START_REF] Isik | Genomic selection in forest tree breeding: the concept and an outlook to the future[END_REF][START_REF] Van Nocker | Breeding better cultivars, faster: applications of new technologies for the rapid deployment of superior horticultural tree crops[END_REF]. In oil palm, previous studies showed that GS could allow selecting individuals without progeny tests [START_REF] Wong | Genomewide selection in oil palm: increasing selection gain per unit time and cost with small populations[END_REF]Cros et al. 2015b). However, GS also has the potential to improve the analysis of progeny tests. So far, no empirical study has investigated whether GS can increase the accuracy of the GCA of progeny-tested oil palms.

The GS model G-BLUP [START_REF] Vanraden | Genomic measures of relationship and inbreeding[END_REF][START_REF] Habier | The impact of genetic relationship information on genome-assisted breeding values[END_REF]) is a mixed model that makes use of molecular information through a genomic matrix (G) of realized relationships. Multivariate analysis using mixed modeling (i.e. considering several traits jointly) aims to take advantage of the genetic correlation between traits to increase accuracy [START_REF] Gilmour | ASReml user guide release 3[END_REF]. Simulations have shown that multivariate G-BLUP can yield higher accuracy than univariate G-BLUP, depending on the heritability of the traits (h²) and their genetic correlation [START_REF] Calus | Accuracy of multi-trait genomic selection using different methods[END_REF][START_REF] Jia | Multiple Trait Genomic Selection Methods Increase Genetic Value Prediction Accuracy[END_REF][START_REF] Guo | Comparison of single-trait and multiple-trait genomic prediction models[END_REF]. When considering two traits with different h², bivariate G-BLUP led to a higher increase in accuracy for the trait with the lowest h². In addition, [START_REF] Jia | Multiple Trait Genomic Selection Methods Increase Genetic Value Prediction Accuracy[END_REF] and [START_REF] Calus | Accuracy of multi-trait genomic selection using different methods[END_REF] found that the stronger the genetic correlation, the greater the benefit of using a bivariate G-BLUP. However, [START_REF] Jia | Multiple Trait Genomic Selection Methods Increase Genetic Value Prediction Accuracy[END_REF] did not provide evidence for improved accuracy with multivariate analysis when they used empirical data. In oil palm simulations, Cros et al (2015a) used multivariate models but did not make comparisons with univariate models.

The aim of the present study was to compare the potential of univariate and bivariate G-BLUP with that of the current univariate pedigree based BLUP for the analysis of BN and ABW traits in oil palm progeny tests using real data.

Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016) h² ABW = 0.23 ± 0.04) and lower than h² ABW in B (h² BN = 0.5 ± 0.05, h² ABW = 0.57 ± 0.04) (Cros et al. 2015b).

Molecular data

Among the progeny-tested parents, 140 group A and 131 group B individuals were genotyped.

Supplementary Table S1 lists the distribution of these individuals among the populations constituting the parental groups. Genotyping was performed with 313 simple sequence repeat markers (SSR) [START_REF] Billotte | Microsatellite-based high density linkage map in oil palm (Elaeis guineensis Jacq.)[END_REF][START_REF] Tranbarger | SSR markers in transcripts of genes linked to post-transcriptional and transcriptional regulatory functions during vegetative and reproductive development of Elaeis guineensis[END_REF][START_REF] Zaki | Elaeis oleifera Genomic-SSR Markers: Exploitation in Oil Palm Germplasm Diversity and Cross-Amplification in Arecaceae[END_REF]. Phenotypic observation of the fruit type (i.e. shell thickness) was included as a two-allele marker, corresponding to the Sh gene (Singh et al 2013). Missing data (1.7% in group A and 2.9% in group B) were imputed with BEAGLE 3.3.2 [START_REF] Browning | Rapid and Accurate Haplotype Phasing and Missing-Data Inference for Whole-Genome Association Studies By Use of Localized Haplotype Clustering[END_REF]. Finally, group A had 265 polymorphic SSR (mean 3.05 alleles ± 0.89 (standard deviation)), and group B had 289 polymorphic SSR (mean 6.25 alleles ± 2.35). For each group only the polymorphic markers were used for the genomic models.

Prediction models

Univariate T-BLUP

The traditional pedigree-based mixed model or T-BLUP was used to predict the genetic effects, i.e. the general combining abilities (GCAs) in A x B crosses of the progeny-tested parents and the specific combining abilities (SCAs) of the crosses (dominance effects). The model was:

P = Xβ + Wb + Z A g A + Z B g B + Z D s + e [1]
where P is the vector of hybrid phenotypes (BN or ABW), X and W are incidence matrices of the experimental design effects, β and b are the vectors of fixed and random effects due to the experimental design, respectively,

Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016) The random genetic effects followed the model of [START_REF] Stuber | Gene effects and variances in hybrid populations[END_REF], with g A ~ N(0, [START_REF] Stuber | Gene effects and variances in hybrid populations[END_REF][START_REF] Lynch | Genetics and analysis of quantitative traits[END_REF].

𝜎𝜎 𝑔𝑔 A 2 × A A ), g B ~ N(0, 𝜎𝜎 𝑔𝑔 B 2 × A B )
Fixed effects were: overall mean, "trial" (26 levels), "block" (152 levels) and "age" (6 levels). Random effects associated with the experimental design were "elementary plots" (3,464 levels), "individual" (30,872 levels), interaction "age*cross" ("α*s", 2,855 levels); with "individual" nested in "elementary plots", "elementary plots" nested in "block", and "block" nested in "trial". The random experimental design effects followed a normal distribution of the form N(0, σ² × I), where I is the identity matrix and σ² the associated variance, with the exception of α*s that followed N(0, 𝜎𝜎 𝛼𝛼 * 𝑠𝑠 2 × I 6×6 ⊗ D). The errors e followed N(0, σ² e × I 180872×180872 ), where σ² e is the residual variance.

Variance parameters were estimated by restricted maximum likelihood (REML) and solutions of the mixed model were obtained by resolving Henderson's mixed model equations [START_REF] Henderson | Best Linear Unbiased Estimation and Prediction under a Selection Model[END_REF], using R-ASReml [START_REF] Gilmour | ASReml user guide release 3[END_REF][START_REF] Core | R: a language and environment for statistical computing[END_REF].

Univariate G-BLUP

In the genomic selection model G-BLUP, the pedigree coancestry matrices used in [1] were replaced by additive genomic coancestry matrices G A and G B for groups A and B, respectively.

As some progeny-tested individuals were not genotyped, their pedigree coancestry had to be combined with the molecular coancestry of the genotyped individuals. For this purpose, we used the single-step approach with matrices H A , H B and D H designed to combine both genomic and pedigree information [START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF][START_REF] Christensen | Genomic prediction when some animals are not genotyped[END_REF]. For each parental group, H inverse was built as follows:

𝐇𝐇 -𝟏𝟏 = 𝐀𝐀 -𝟏𝟏 + � 𝟎𝟎 𝟎𝟎 𝟎𝟎 𝐆𝐆 -𝟏𝟏 -𝐀𝐀 𝟐𝟐𝟐𝟐 -𝟏𝟏 �
Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016) [START_REF] Lynch | Estimation of relatedness by DNA fingerprinting[END_REF][START_REF] Li | Similarity of DNA fingerprints due to chance and relatedness[END_REF]) and was defined as:

𝐆𝐆 𝐀𝐀𝐀𝐀𝐀𝐀 = 𝐙𝐙𝐙𝐙 𝐭𝐭 4L
where Z is the genotypic matrix with as many columns as alleles, with the individuals in rows, and containing in the i th column the number of copies of the i th allele (Z xy ⋲ {0,1,2}), and L is the total number of markers. This index estimates coancestry from alike-in-state (AIS) alleles, and assumes that each allele was unique in the founder population that generated the population under study [START_REF] Eding | Marker-based estimates of between and within population kinships for the conservation of genetic diversity[END_REF].

G OF was obtained according to [START_REF] Vanraden | Genomic measures of relationship and inbreeding[END_REF][START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF], with a modification for multiallelic markers:

𝐆𝐆 𝐎𝐎𝐎𝐎 = (𝐙𝐙 -𝐏𝐏)(𝐙𝐙 -𝐏𝐏) t 4 ∑ (1 -∑ ∑ 𝑝𝑝 𝑙𝑙 𝑎𝑎 2 𝑎𝑎 𝑙𝑙 L l=1
) where P is a matrix containing twice the observed allelic frequency (OF) of the i th allele in the genotyped individuals in the i th column.

The coancestry matrix of [START_REF] Vanraden | Genomic measures of relationship and inbreeding[END_REF][START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF] normally requires the allele frequencies in the founder population. As these frequencies are usually not known, they are commonly replaced by the observed frequencies, as we did in our study. G N was derived from G OF , with normalization to provide more realistic variance and accuracy estimations [START_REF] Forni | Different genomic relationship matrices for single-step analysis using phenotypic, pedigree and genomic information[END_REF]:

𝐆𝐆 𝐍𝐍 = 1 2 × (𝐙𝐙 -𝐏𝐏)(𝐙𝐙 -𝐏𝐏) t {trace[(𝐙𝐙 -𝐏𝐏)(𝐙𝐙 -𝐏𝐏) t ]} n ⁄
where n is the number of genotyped individuals.

Multivariate models

Multivariate models were built from [1], as follows (Mrode 2005):

Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016)

� 𝐀𝐀𝐀𝐀𝐀𝐀 𝐀𝐀𝐍𝐍 � = � 𝐗𝐗 𝟎𝟎 𝟎𝟎 𝐗𝐗 � � 𝛃𝛃 ABW 𝛃𝛃 BN � + � 𝐀𝐀 𝟎𝟎 𝟎𝟎 𝐀𝐀 � � 𝐛𝐛 ABW 𝐛𝐛 BN � + � 𝐙𝐙 A 𝟎𝟎 𝟎𝟎 𝐙𝐙 A � � 𝐠𝐠 𝐀𝐀 𝐀𝐀𝐀𝐀𝐀𝐀 𝐠𝐠 𝐀𝐀 𝐀𝐀𝐍𝐍 � + � 𝐙𝐙 B 𝟎𝟎 𝟎𝟎 𝐙𝐙 B � � 𝐠𝐠 𝐀𝐀 𝐀𝐀𝐀𝐀𝐀𝐀 𝐠𝐠 𝐀𝐀 𝐀𝐀𝐍𝐍 � + � 𝐙𝐙 D 𝟎𝟎 𝟎𝟎 𝐙𝐙 D � � 𝐬𝐬 𝐀𝐀𝐀𝐀𝐀𝐀 𝐬𝐬 𝐀𝐀𝐍𝐍 � + � 𝐞𝐞 𝐀𝐀𝐀𝐀𝐀𝐀 𝐞𝐞 𝐀𝐀𝐍𝐍 �
In multivariate T-BLUP, genetic effects were structured as:

� 𝐠𝐠 𝐀𝐀 𝐀𝐀𝐀𝐀𝐀𝐀 𝐠𝐠 𝐀𝐀 𝐀𝐀𝐍𝐍 � ~𝑁𝑁(0, � 𝜎𝜎² 𝑔𝑔 A ABW C g A C g A 𝜎𝜎² 𝑔𝑔 A BN � ⊗ 𝐀𝐀 𝐀𝐀 ) � 𝐠𝐠 𝐀𝐀 𝐀𝐀𝐀𝐀𝐀𝐀 𝐠𝐠 𝐀𝐀 𝐀𝐀𝐍𝐍 � ~𝑁𝑁(0, � 𝜎𝜎² 𝑔𝑔 B ABW C g B C g B 𝜎𝜎² 𝑔𝑔 B BN � ⊗ 𝐀𝐀 𝐀𝐀 ) � 𝐬𝐬 𝐀𝐀𝐀𝐀𝐀𝐀 𝐬𝐬 𝐀𝐀𝐍𝐍 � ~𝑁𝑁(0, � 𝜎𝜎² 𝑠𝑠 ABW C s C s 𝜎𝜎² 𝑠𝑠 BN � ⊗ 𝐃𝐃)
where C gA and C gB are additive genetic covariances and C s is the dominance genetic covariance. Residual effects were structured as:

� 𝐞𝐞 𝐀𝐀𝐀𝐀𝐀𝐀 𝐞𝐞 𝐀𝐀𝐍𝐍 � ~𝑁𝑁(0, � 𝜎𝜎² 𝑒𝑒 ABW C e C e 𝜎𝜎² 𝑒𝑒 BN � ⊗ 𝐀𝐀)
For multivariate G-BLUP, A A , A B and D were replaced by H A , H B and D H , respectively.

Non-genetic random effects had unstructured variances-covariances.

Variances and covariances of both non-genetic effects and genetic effects were estimated by REML.

Comparison of models

For a given type of model (i.e. univariate for ABW, univariate for BN and bivariate), the G-BLUP approaches based on the three additive genomic matrices were compared between themselves and with the T-BLUP model. At this stage, for computational reasons, the SCA effects were considered uncorrelated between traits in the multivariate models. Only the additive genetic variance-covariance structure matrix varied, while the number of observations and estimated parameters remained constant. The models were consequently directly compared based on their deviance (-2LogLikelihood), which was the equivalent of comparing their Akaike information criterion and Bayesian information criterion. The convergence of REML algorithm was also considered.

The univariate G-BLUP and multivariate G-BLUP models were compared based on their accuracy, which is the correlation between the predicted genetic effects (GCAs or SCAs) and their true value (unknown).

Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016) The accuracy of the genetic effect predicted for the x th level (i.e. parent for GCA or cross for SCA) was estimated from its relation with the prediction error variance (PEV) [START_REF] Clark | The importance of information on relatives for the prediction of genomic breeding values and the implications for the makeup of reference data sets in livestock breeding schemes[END_REF]:

r 𝑥𝑥 = �1 - PEV 𝑥𝑥 𝜎𝜎 2 𝚺𝚺 𝑥𝑥𝑥𝑥
where σ² is the variance of the genetic effect, Σ xx is the x th term of the diagonal of the associated variancecovariance matrix and PEV 𝑥𝑥 = (u � -u) 𝑥𝑥 2 , with u the genetic effect considered. PEVs were computed from the elements of the inverse of the mixed model equations, based on theoretical derivations from [START_REF] Henderson | Best Linear Unbiased Estimation and Prediction under a Selection Model[END_REF] (i.e. not obtained by cross-validation). Consequently, for any progeny-tested individual x, the accuracy associated with its GCA for a given trait was:

r GCA 𝑥𝑥 = � 1 - PEV GCA 𝑥𝑥 𝐆𝐆 𝑥𝑥𝑥𝑥 𝜎𝜎 𝑔𝑔 2 [3]
where 𝜎𝜎 𝑔𝑔 2 is the estimated additive variance of the trait for the parental group of x. For the univariate and multivariate G-BLUP and for each trait, we computed the mean r GCA over the 140 group A parents and the 131 group B parents that were genotyped. For each group and each trait, the mean r GCA of univariate and multivariate models was compared using a t-test and a Bonferroni correction. We also compared the accuracy of SCA, which for cross x × y and a given trait, was:

r SCA 𝑥𝑥𝑥𝑥 = �1 - PEV SCA 𝑥𝑥𝑥𝑥 𝐃𝐃 𝑥𝑥𝑥𝑥 𝑥𝑥𝑥𝑥 σ 𝑠𝑠 2 [4]
where 𝐃𝐃 𝑥𝑥𝑥𝑥 𝑥𝑥𝑥𝑥 = G 𝑥𝑥𝑥𝑥 G 𝑥𝑥𝑥𝑥 and 𝜎𝜎 𝑠𝑠 2 the estimated dominance variance for the trait. For the univariate and multivariate G-BLUP and for each trait, we computed the mean r SCA over the 478 crosses evaluated in the progeny test and over 256 crosses that had not been evaluated. These 256 crosses were sampled from the unevaluated crosses among all possible crosses between the 140 A and 131 B parents, with a balanced representation of the parents of both groups (i.e. each parent occurred once or twice among the 256 unevaluated crosses). To obtain the PEV SCA of the 256 unevaluated crosses and compute their r SCA , these crosses were added to the 𝐃𝐃 matrix prior to analyzing the mixed models, following [START_REF] Henderson | Best linear unbiased prediction of breeding values not in the model for records[END_REF]. The mean r SCA of the univariate and multivariate models was compared using a t-test for each group and trait, and a Bonferroni correction was applied to adjust the p-values.

Marker density

Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016) We studied the effect of marker density on the prediction of GCAs by the multivariate G-BLUP model.

This was investigated independently in the two parental groups by varying the number of markers for one group, while keeping the maximum number of markers for the other group. The number of markers m varied from 10 to 265 in group A and from 10 to 289 in group B, with a step of 10. At each density, five replicates were made, for each replicate, we used a random subset of m markers chosen among all the available polymorphic markers for the group. For each replicate, the additive coancestry matrix of the group concerned was calculated using the m markers, and the dominance matrix was calculated using the m markers for the group concerned and all the markers of the other group. To assess the effect of the number of markers on the prediction of GCAs, we calculated the prediction accuracy of the model, i.e. the correlation between the predicted GCAs (for the group whose marker density varied) and the reference GCAs. The reference GCAs were obtained from the most accurate model previously identified (actually the multivariate G-BLUP) using all the markers, so that the prediction accuracy was the best approximation of accuracy. The different levels of marker number, of replicates per level of marker number and the two parental groups meant the calculations had to be repeated many times, so, to speed up the process, no covariance was specified for the dominance effects in the multivariate G-BLUP model used here.

Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016)

Coancestry matrices

The distribution of coancestry estimates in group A and group B is shown in Figure 1. Coancestry estimates in G AIS and A belonged to [START_REF] Grattapaglia | Breeding forest trees by genomic selection: current progress and the way forward[END_REF]1], as expected, as coancestry is the probability that two alleles on a random locus of two individuals are identical by descent [START_REF] Wright | Coefficients of inbreeding and relationship[END_REF][START_REF] Malécot | Les mathématiques de l'hérédité[END_REF]. The median value of the two VanRaden matrices (G OF , G N ) was below 0, meaning that more than half the coancestry estimates were negative. The REML algorithm converged with A and G AIS matrices. The smallest deviance was obtained with G AIS (Table 1). The G OF and G N matrices were not positive definite and the REML algorithm did not converge, leading to higher deviances than with G AIS and A. Therefore, for our dataset, the G AIS matrix appeared to be more appropriate than the other genomic matrices G OF and G N , and than the genealogical matrix A. For the rest of the study, we consequently only used G AIS in the G-BLUP.

For both A and B groups, coancestry estimates in G AIS were higher than in matrix A. G AIS did not contain any null coancestry estimates, whereas A contained 73.6% null coancestry estimates for group A and 42.9% for group B. The coancestry estimates for group A were lower than those for group B in the A matrix, but were higher in G AIS . The variability in coancestry estimates was higher in group A than in group B.

Multivariate G-BLUP

The multivariate G-BLUP revealed very high additive correlations, reaching -0.997 in the parental group A and -0.917 in group B, very high dominance correlations (-0.987) and low residual correlations (-0.158).

The GCA accuracy of the univariate and multivariate G-BLUP are depicted in Figure 2A. For all combinations of groups and traits, mean GCA accuracy was higher with the multivariate G-BLUP model than with the univariate G-BLUP (p<10 -100 ). According to the Pearson correlation, the GCAs obtained with the multivariate T-BLUP and univariate G-BLUP models were highly correlated with the reference GCAs (average Pearson correlation coefficient of 0.978 and 0.966, respectively). However, the Spearman's correlation coefficients computed on the top 10% individuals

were not as high, with an average value of 0.595, ranging from 0.213 to 0.978. This indicated that the model impacted the selection of the progeny tested individuals, and was therefore of importance for practical breeding.

In addition, the multivariate T-BLUP gave GCAs with ranks that were more correlated with the ranks of the reference GCAs than the univariate G-BLUP (Spearman's rank correlation coefficient of 0.696 and 0.562, respectively). Therefore, the improvement obtained in the GCA estimates when using a multivariate genomic approach compared to the conventional univariate T-BLUP resulted more from the multivariate analyze than from the use of the genomic data.

SCA accuracy was higher with the multivariate G-BLUP than with the univariate G-BLUP (p<10 -100 ) (Figure 2B andC). The average increase in SCA accuracy was 18.7%, ranging from 12.9% (trait BN, unevaluated crosses) to 24.6% (trait ABW, evaluated crosses). With the multivariate G-BLUP model, SCA accuracies were on average 0.76 for evaluated crosses and 0.68 for unevaluated crosses.

The h² obtained with the multivariate genomic model were h² BN = 0.53 and h² ABW = 0.35 in group A, and h² BN = 0.4 and h² ABW = 0.79 in group B (see Supplementary Table S2 for the detail of variances).

Marker density

Figure 3 shows the effect of marker density on the prediction accuracy of GCAs with the multivariate G-BLUP for ABW. The results obtained for BN were very similar (Supplementary Fig. S1), certainly due to the high genetic correlation between ABW and BN. As marker density increased, the prediction accuracy of multivariate G-BLUP also rapidly increased before reaching a plateau slightly above the prediction accuracy of multivariate T-BLUP. To outperform the prediction accuracy of the multivariate T-BLUP model, multivariate G-BLUP required 110 markers for group A and 70 markers for group B, for both ABW and BN traits. The prediction accuracies exceeded 0.99 with 160 markers for group A for both ABW and BN, and with 80 (respectively 90) markers for group B for ABW (respectively BN).

Discussion

Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016) The general combining ability (GCA) for bunch number (BN) and average bunch weight (ABW) of progeny-tested oil palms is currently obtained with a pedigree-based univariate mixed model analysis of phenotypic data of hybrid individuals. In this study, we showed that using a multivariate model and replacing the genealogical coancestry matrices by molecular matrices of realized coancestry improved the analysis, leading to better estimated GCAs. In addition, the accuracy of the SCAs, usually neglected, reached interesting levels. We also found that this could be achieved with a reduced marker density. Indeed, the number of SSR markers that enabled G-BLUP to reach the same prediction accuracy as T-BLUP was 110 for group A and 70 for group B;

while 160 markers in group A and 90 in group B were needed to achieve the maximum benefit offered by the genomic approach.

Genomic versus genealogical coancestries

We observed many null coancestry estimates in A, whereas all coancestry estimates in G were higher than zero. This reflected the fact that the pedigrees used to estimate the A matrices did not reach the base of the unrelated founders of the different populations. Consequently, the pedigree-based coancestries underestimated the real coancestries, whereas the genomic coancestries were able to capture these hidden relationships, which did not appear in the pedigree. However, as G AIS considered identity by state and A identity by descent, the values in G AIS were actually overestimated if several copies of some alleles were present in the founder populations [START_REF] Eding | Marker-based estimates of between and within population kinships for the conservation of genetic diversity[END_REF]. Nevertheless, the G-BLUP model using G AIS was more appropriate for the data than the T-BLUP model, as shown by its higher likelihood.

The Deli individuals, which made up most of group A, originated from four oil palms planted in 1848 in Indonesia, while the African populations in group B can be traced back to the first half of the 20 th century, with around 15 to 20 founders [START_REF] Corley | Selection and breeding[END_REF]. The higher G AIS values found in group A than in group B is not surprising, given the longer history of inbreeding, drift and artificial selection of Deli individuals. However, the coancestry estimates for group A were lower than those for group B in the A matrix. This resulted from the depth (number of generations) of the pedigree and from the history of the populations constituting the parental groups. In group B, the data available on the pedigrees referred roughly to the initial generation, but the longer history of the Deli population was not covered by its pedigree, which did not go back far enough in time. This explained why, according to the pedigrees, there were fewer coancestries in group A than in group B. This also explained the fact that the number of relationships hidden in the pedigree but captured by the markers was higher in group A than in group B. This increased the benefit of using the genomic models more for group A than for

Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016)

Multivariate model

This is the first study to investigate the benefit of using multivariate genomic models for oil palm breeding. Using empirical data, we demonstrated that multivariate genomic models improved the prediction accuracy of additive effects (GCAs). In addition, we showed that in each parental group, the trait with the lowest heritability (ABW in group A and BN in group B) benefited the most from the use of a multivariate model. Both findings are in agreement with the results of previous simulations [START_REF] Calus | Accuracy of multi-trait genomic selection using different methods[END_REF][START_REF] Jia | Multiple Trait Genomic Selection Methods Increase Genetic Value Prediction Accuracy[END_REF][START_REF] Guo | Comparison of single-trait and multiple-trait genomic prediction models[END_REF] but, in addition to the results of these studies, we showed that genomic multivariate models also increased the prediction accuracy of dominance effects (SCAs).

In the multivariate G-BLUP model, covariance between traits is considered to be identical at each marker. This could reduce the efficiency of multivariate G-BLUP relatively to a multivariate Bayesian method that would allow marker specific covariances between traits [START_REF] Guo | Comparison of single-trait and multiple-trait genomic prediction models[END_REF]). An empirical comparison of these two statistical approaches with oil palm data would thus be useful.

Density and type of molecular markers

In the conventional pedigree-based analysis of progeny tests, the GCA of a progeny-tested individual results from the phenotypes of its progeny and the progeny of its relatives. The measure of coancestry used in this conventional approach is an expected value, as it is based on pedigree, and may thus differ from the true coancestry. The genomic approach improves this situation as it uses the realized coancestry between progenytested individuals. We found that even small numbers of markers (110 in group A and 70 in group B) gave

GCAs similar to those obtained with a conventional pedigree-based model, which was likely a consequence of the small effective size of the parental groups of oil palm (<10) (Cros et al. 2015b). The respective history of the parental groups, with the longer history of inbreeding, drift and artificial selection in group A than group B, led to less variable realized coancestries in group A and, as we used SSR markers instead of SNPs as in the other GS studies, to fewer alleles per marker in group A. As a consequence, group A required more markers to capture the realized coancestries than group B. This difference between groups therefore resulted from their contrasted history and from the use of a multiallelic type of markers.

Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016) When all the markers were used, we achieved higher prediction accuracy than with the pedigree-based model. However, the improvement was very limited thus indicating that the phenotypic data of the hybrid progenies play a major role in the quality of the estimation of the GCAs, while the coancestry matrices used in the model play a secondary role.

In the present study, the progeny-tested individuals were genotyped using SSR markers. This type of marker is suitable for genotyping a relatively small number of individuals with a low coverage of the genome, but the practical application of GS for breeding implies large scale genotyping capabilities, at reasonable cost.

Therefore, future GS studies in oil palm will likely use SNP markers, as this would reduce the cost per data point and speed up the genotyping process. Although more SNPs are needed to achieve the same GS accuracy as with SSR markers [START_REF] Solberg | Genomic selection using different marker types and densities[END_REF], the efficiency of the current genotyping technologies ensures that SNP density will not be a limiting factor to implement GS in oil palm. Thus, two SNP arrays have been developped

for this species, with 4.5 K [START_REF] Ting | High density SNP and SSR-based genetic maps of two independent oil palm hybrids[END_REF]) and 200K SNP [START_REF] Teh | Genome-wide association study of oil palm mesocarp oil yield content and its application for marker selection[END_REF]; while [START_REF] Pootakham | Genome-wide SNP discovery and identification of QTL associated with agronomic traits in oil palm using genotypingby-sequencing (GBS)[END_REF] identified over 20 K SNP using the genotyping-by-sequencing approach.

Comparison of models

We must emphasize that, as shown by formulas [3] and [4], accuracy based on prediction error variances (PEV) cannot be used to compare models that differ in their genetic variance covariance matrices, as the estimated variances refer to a different base population. Here, the base population implicitly used with the T-BLUP model was made up of the individuals with no known parents in the pedigrees, while with the G-BLUP model, it was made up of genotyped individuals. In other words, even for methods that yield the same estimated breeding values, accuracies obtained from the PEVs are not invariant to parameterization [START_REF] Stranden | Allele coding in genomic evaluation[END_REF]. Consequently, the fact that the accuracies we obtained from the PEVs for the G-BLUP models were lower than the accuracies of T-BLUP, which were around 0.90 (Cros et al. 2015b), was not meaningful. When evaluating the potential of GS to predict the GCA of individuals that have not been progenytested, accuracy is often estimated using a cross-validation approach, like that used in Cros et al. (2015b) for oil palm. However, in the present study, this was not possible as we were interested in the ability of GS to predict the GCA of progeny-tested individuals, and so the approach we chose was to compare T-BLUP and G-BLUP models based on their likelihood, and to trust the best model. Similarly, when considering either G-BLUP or T-BLUP, it was not possible to use likelihood to compare the univariate and multivariate versions of the model, as

Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016)

Implications for breeding

The choice of the model to analyze the progeny tests impacted the practical breeding work, as it affected the ranking of the evaluated individuals and therefore the set of the selected individuals. Here, we focused on BN and ABW, two major traits determining oil yield, but genomic models could be used instead of the traditional pedigree based models for all the traits recorded in progeny tests, i.e. bunch quality, height increment, disease symptoms [START_REF] Corley | Selection and breeding[END_REF][START_REF] Durand-Gasselin | Sélection du palmier à huile pour une huile de palme durable et responsabilité sociale[END_REF], annual profile of bunch production measured by the Gini coefficient [START_REF] Cros | Differential response to water balance and bunch load generates diversity of bunch production profiles among oil palm crosses (Elaeis guineensis)[END_REF], etc. In addition, correlated traits should be analyzed jointly in a multivariate model. Here we considered two traits but a higher number of correlated traits could easily be used.

In oil palm, several traits are known to be correlated including the number of fruits per bunch and the average fruit weight, the percentage of pulp and the percentage of kernel in the fruits. As indicated by the literature, the benefit of using a multivariate approach will result from the h² of the traits included in the model and from their correlation. Using the same dataset, Cros et al. (2015b) showed that GS could predict the GCA of non-progenytested individuals for some traits in group B, in particular when the candidate individuals were highly related to the training set. Here, we showed that GS was also useful to predict the GCA of progeny-tested individuals and the SCA of crosses. GS is therefore a highly valuable method for oil palm breeding, even with low marker density.

Our experimental design involved a mean number of 65 hybrid individuals per cross. It would be interesting to study the effect of decreasing the number of hybrid individuals per cross in the progeny tests, as we would expect the G-BLUP model to be less affected than the T-BLUP, thanks to the extra information provided to the G-BLUP (realized coancestries). Reducing the number of hybrid individuals per cross would also allow progeny-testing more parents, thus increasing the selection intensity without increasing the cost of the progeny tests. The importance of hybrid phenotypes in the prediction of GCAs also suggests that the number of markers required to predict the GCAs of non-progeny-tested individuals might be higher than the number required for progeny-tested individuals. However, this point requires further investigation.

To our knowledge, this is the first report of accuracy of SCA for oil palm crosses. It appeared to be lower than the accuracy of the GCAs, with the mean accuracy of SCA of crosses that were not evaluated in the fields reaching 0.68 with the multivariate model. The low proportion of dominance variance in total genetic

Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016) Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016)

  and s ~ N(0, 𝜎𝜎 𝑠𝑠 2 × D), where σ² gA and σ² gB are the additive variances of the A and B parents in A × B hybrid crosses, respectively, and σ² s is the variance of the dominance effects in the A × B population. Given the hybrid nature of the crosses, the A matrices contain Malécot's coefficient of coancestry f (Malécot 1948), such as A xy = { f xy } between individuals x and y. They were built from the pedigrees with the R package synbreed (Wimmer et al. 2012). The D matrix is the dominance coancestry matrix between crosses, obtained as D = A A ⊗ A B [2], i.e. with elements D AB,A'B' = f AA' f BB' , as A and B individuals are unrelated

Figure legends Figure 1 Figure 2 Figure 3 Figure 1 Figure 2 Figure 3

 legends123123 Figure legends

Material and methods Experimental population and phenotypes

  Version définitive du manuscrit publié dans / Final version of the manuscript published in :The progeny test involved 146 group A parents crossed with 155 group B parents according to an incomplete factorial design with 478 crosses. The crosses were planted between 1995 and 2000 in 26 trials located in the same area, at the Aek Loba estate (SOCFINDO, North Sumatra). All the vegetal material belonged to the PalmElit breeding program (www.palmelit.com). Annual bunch production data, i.e. bunch number (BN) and average bunch weight (ABW), were collected on 30,872 progeny palms of type tenera (thin-shelled commercial type) from 6 years old up to 11 years old. More details on the experimental design are given inCros et al. (2015b). Phenotypic correlation between ABW and BN was -0.682. The narrow-sense heritabilities h² of ABW and BN varied with the parental population, h² BN was higher than h² ABW in A (h² BN = 0.31 ± 0.04 [s.e.],

		Molecular Breeding, 2016, 36(1), 36:2	http://dx.doi.org/10.1007/s11032-015-0423-1
	Manuscrit d'auteur / Author Manuscript	2.
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript

. Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1

  Version définitive du manuscrit publié dans / Final version of the manuscript published in : A , Z B and Z D are incidence matrices of the genetic random effects, g A and g B are the vectors of GCA of parents A and B, respectively, s is the vector of SCA of crosses, and e is the vector of residual effects.

	Molecular Breeding, 2016, 36(1), 36:2	http://dx.doi.org/10.1007/s11032-015-0423-1
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript

. Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1 Z

AIS , G OF and G N . G AIS used a similarity index

		Version définitive du manuscrit publié dans / Final version of the manuscript published in :
		Molecular Breeding, 2016, 36(1), 36:2	http://dx.doi.org/10.1007/s11032-015-0423-1
	Manuscrit d'auteur / Author Manuscript	where A is the pedigree coancestry matrix including all the individuals in the group, A 22 and G are the pedigree and the genomic coancestry matrices, respectively, containing only the genotyped individuals. Then H D was built in the same way as in equation [2]: D H = H A ⊗ H B This led to g A ~ N(0, σ² gA × H A ), g B ~ N(0, σ² gB × H B ) and s ~ N(0, σ² s × D H ). Three different genomic additive coancestry matrices G were compared: G
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript
			. Multivariate genomic
		model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular
		Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1

  Version définitive du manuscrit publié dans / Final version of the manuscript published in :

	Molecular Breeding, 2016, 36(1), 36:2	http://dx.doi.org/10.1007/s11032-015-0423-1
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript

. Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1

  Version définitive du manuscrit publié dans / Final version of the manuscript published in :

		Molecular Breeding, 2016, 36(1), 36:2	http://dx.doi.org/10.1007/s11032-015-0423-1
	Manuscrit d'auteur / Author Manuscript	3. Results
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript

. Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1

  The average increase was 22.5%, ranging from 13.2% for ABW in group B to 32.1% for BN in group B. There were differences in GCA accuracy between traits within a parental group with univariate G-BLUP, but the multivariate G-BLUP model increased the GCA accuracy of both traits to the same level, i.e. 0.83 in group B and 0.88 in group A. Thus, the trait with the lowest GCA accuracy in the univariate models (ABW for group A and BN for group B) benefited the most from the multivariate model.As the multivariate G-BLUP model predicted GCAs best, we used the GCAs predicted by this model as reference GCAs. The Pearson correlation coefficients between GCAs predicted by any of the models (univariate Version définitive du manuscrit publié dans / Final version of the manuscript published in : ). The GCAs obtained with univariate T-BLUP were generally the least correlated with the reference GCAs, with an average Pearson correlation coefficient of 0.946 and Spearman's correlation coefficient of 0.527.

	Molecular Breeding, 2016, 36(1), 36:2	http://dx.doi.org/10.1007/s11032-015-0423-1
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript

or multivariate, T-BLUP or G-BLUP) and reference GCAs are listed for each trait and each group in Table

2

, as well as the Spearman's rank correlation of the 10% best individuals ("best" when evaluated by the reference Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M.,

Cros, D. (2016)

. Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1 model

  Version définitive du manuscrit publié dans / Final version of the manuscript published in :

	Molecular Breeding, 2016, 36(1), 36:2	http://dx.doi.org/10.1007/s11032-015-0423-1
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript

. Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1

  Version définitive du manuscrit publié dans / Final version of the manuscript published in : B, as shown by the bigger increase in the correlation with reference GCAs in group A when the G-BLUP model was used instead of the T-BLUP, than in group B.

	Molecular Breeding, 2016, 36(1), 36:2	http://dx.doi.org/10.1007/s11032-015-0423-1
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript

. Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1 group

  Version définitive du manuscrit publié dans / Final version of the manuscript published in :

	Molecular Breeding, 2016, 36(1), 36:2	http://dx.doi.org/10.1007/s11032-015-0423-1
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript

. Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1

  Version définitive du manuscrit publié dans / Final version of the manuscript published in : datasets (phenotypic observations) differed, but using PEV-based accuracy was relevant as the variance covariance matrices were the same for the univariate and multivariate models.

	Molecular Breeding, 2016, 36(1), 36:2	http://dx.doi.org/10.1007/s11032-015-0423-1
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript

. Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1 the

Table 1

 1 Version définitive du manuscrit publié dans / Final version of the manuscript published in : Deviance of the mixed model according to the coancestry matrices (G AIS ¸ G OF , G N and A), for average

	Molecular Breeding, 2016, 36(1), 36:2	http://dx.doi.org/10.1007/s11032-015-0423-1

. Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1

Table 2

 2 Pearson correlation (top) and Spearman's rank correlation of the top 10% individuals (bottom) between predicted GCAs produced by a G-BLUP or T-BLUP, univariate or multivariate model and the reference GCAs from the multivariate G-BLUP. The Pearson correlation coefficients were calculated based on the 140 group A genotyped parents and the 131 group B genotyped parents. The top 10% individuals represented 14 individuals Version définitive du manuscrit publié dans / Final version of the manuscript published in : Version définitive du manuscrit publié dans / Final version of the manuscript published in :

	in group A and 13 in group B
	Pearson correlation		Group A	Group B
			ABW	BN	ABW	BN	mean
	Multivariate	G-BLUP T-BLUP	1 0.971	1 0.971	1 0.986	1 0.983 0.978 1
	Univariate	G-BLUP T-BLUP	0.963 0.905	0.980 0.960	0.976 0.969	0.946 0.966 0.947 0.946
	Spearman's rank correlation on the top 10% individuals	Group A	Group B
			ABW	BN	ABW	BN	mean
	Multivariate	G-BLUP T-BLUP	1 0.732	1 0.424	1 0.978	1 0.648 0.696 1
	Univariate	G-BLUP T-BLUP	0.789 0.635	0.218 0.213	0.830 0.890	0.412 0.562 0.368 0.527

Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M.,

Cros, D. (2016)

. Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1 Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016). Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1

  Version définitive du manuscrit publié dans / Final version of the manuscript published in : Molecular Breeding, 2016, 36(1), 36:2 http://dx.doi.org/10.1007/s11032-015-0423-1

	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript

. Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1

Table S1

 S1 Details on the 271 parents used in the study, per group and population. All these individuals were present in the pedigree and genotyped.Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M.Version définitive du manuscrit publié dans / Final version of the manuscript published in : Molecular Breeding, 2016, 36(1), 36:2 http://dx.doi.org/10.1007/s11032-015-0423-1

	Manuscrit d'auteur / Author Manuscript	Group A B	Population Deli Angola Total La Mé Yangambi Nigeria La Mé × Yangambi La Mé × Sibiti Total	Total 131 9 140 93 24 2 5 7 131
	Manuscrit d'auteur / Author Manuscript
	Manuscrit d'auteur / Author Manuscript

, Cros, D. (2016). Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests. Molecular Breeding, 36 (2), Non paginé. DOI : 10.1007/s11032-015-0423-1

Table S2

 S2 Variances estimated with the multivariate genomic model for average bunch weight (ABW) and bunch number (BN): additive variances for parental groups A (𝜎𝜎² 𝑔𝑔 A ) and B (𝜎𝜎² 𝑔𝑔 B ) and dominance variance (𝜎𝜎² 𝑠𝑠 ) in A x B crosses

		𝜎𝜎² 𝑔𝑔 A	𝜎𝜎² 𝑔𝑔 B	𝜎𝜎² 𝑠𝑠
	ABW	1.15	2.62	2.81
	BN	3.11	2.34	7.99

Acknowledgements

We acknowledge SOCFINDO (Indonesia) and CRAPP (Benin) for planning and carrying out the field trials with CIRAD (France) and authorizing the use of the data for this study. We thank P. Sampers and V.

Pomiès for help in genotyping. This research was partly funded by a grant from PalmElit SAS.

Comment citer ce document : Marchal, A., Legarra , A., Tisne, S., Carasco-Lacombe, C., Manez, A., Suryana, E., Omoré, A., Nouy, B., Durand-Gasselin, T., Sanchez, L., Bouvet, J.-M., Cros, D. (2016) [START_REF] Soh | Ranking parents by best linear unbiased prediction (BLUP) breeding values in oil palm[END_REF], many breeding companies have not yet started using BLUP for practical breeding decisions and, those that have started, did so relatively recently. Without BLUP taking coancestries into account, the analysis of the progeny-tests only provides SCA estimates for the crosses that were evaluated. When the BLUP model is provided with pedigree information, most of the dominance relationship matrix D contains zero elements due to the numerous null coancestries in the A matrices. In these conditions, the BLUP model will yield no estimates of SCA at all or only very inaccurate estimates for crosses that were not evaluated in the trials.

However, as markers are more efficient than pedigrees at capturing coancestries, the accuracy of SCAs obtained with GBLUP is high enough to make selection possible, particularly with multivariate analysis of correlated traits. For these reasons, oil palm breeders should also consider SCAs when selecting among progeny-tested parents, since, although relatively small, the extra genetic gain obtained compared to selection based only on GCAs, would come at no extra cost.

Conflicts of interest

The authors declare no conflict of interest.