David Cros
email: david.cros@cirad.fr

Luther Mbo-Nkoulou

Joseph Martin Bell

Jean Oum

Aurélien Masson

Dinh Mouman Soumahoro

Minh Tran

Zeineb Achour

Vincent Le Guen

André Clément-Demange

Joseph Martin

André Clement-Demange

Within-family genomic selection in rubber tree (Hevea brasiliensis) increases genetic gain for rubber production

Keywords: marker assisted selection, genomic predictions, selection response, clonal varieties

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

Introduction

The rubber tree (Hevea brasiliensis, hereafter Hevea) is almost the only source of commercial natural rubber (1,4 cis-polyisoprene), 70% of which is used by the tire industry. The production of natural rubber worldwide has increased steadily over time, and is now exceeding 12 Mt yearly (FAOSTAT, 2018). The total cultivated area, currently over 11 million hectares, is held by smallholders (80%) and industrial estates (20%). More than 90% of the production takes place in Asia, with Thailand and Indonesia as the largest producers. Côte d'Ivoire is the seventh world producer, and produced 420 thousand tons in 2017. Predictions indicate that demand for natural rubber will exceed 19 Mt in 2025 [START_REF] Warren-Thomas | Increasing Demand for Natural Rubber Necessitates a Robust Sustainability Initiative to Mitigate Impacts on Tropical Biodiversity[END_REF], even though rubber plantations are already responsible for deforestation and pose threats to biodiversity, in particular in South-East Asia [START_REF] Ahrends | Current trends of rubber plantation expansion may threaten biodiversity and livelihoods[END_REF][START_REF] Warren-Thomas | Increasing Demand for Natural Rubber Necessitates a Robust Sustainability Initiative to Mitigate Impacts on Tropical Biodiversity[END_REF]. Yield therefore needs to be intensified in existing plantations to meet the expected demand while minimising environmental cost and increasing the income of poor producers.

Genomic selection (GS) [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF], the state-of-the-art method of marker-assisted selection for quantitative traits, can play a key role in taking up this challenge.

Hevea is a diploid species (2n = 36) belonging to the family Euphorbiaceae and originating from the Amazonian forest. Vegetative multiplication by grafting permitted the development of clonal varieties from axillary buds grafted on seedling rootstocks (rubber clones). The initial 'primary' clones derived from ortet selection among populations of non-budded trees resulting from natural pollination. Controlled recombination by hand pollination was then applied to cross the best clonal parents for the generation of full-sib families. However, the naturally-low female fertility of Hevea makes it difficult to construct complex populations of connected families, with highly incomplete mating designs and strong imbalance in family sizes. This generally did not allow to accurately estimate parental genetic values and to take advantage of the large within-families variability. This prompted us to adopt a within-family clonal breeding program which focused specifically on certain large-sized F1 families (≥200 individuals) obtained with the few parent trees that combined good agronomic performances and female fertility sufficient to reach the targeted size. Since the 1990s and the development of molecular genetic markers, these large and highly performing F1 families also gave the opportunity to acquire genetic information about the parents using genetic mapping and quantitative trait loci (QTL) detection [START_REF] Clément-Demange | Hevea Rubber Breeding and Genetics[END_REF]. These strategies of QTL detection have been applied to various traits: resistance to Pseudocercospora ulei [START_REF] Le Guen | A rubber tree's durable resistance to Microcyclus ulei is conferred by a qualitative gene and a major quantitative resistance factor[END_REF], 2007;[START_REF] Lespinasse | Identification of QTLs involved in the resistance to South American leaf blight (Microcyclus ulei) in the rubber tree[END_REF] and to Corynespora cassiicola [START_REF] Tran | Genetic determinism of sensitivity to Corynespora cassiicola exudates in rubber tree (Hevea brasiliensis)[END_REF], vegetative growth and latex production [START_REF] An | A highdensity genetic map and QTL mapping on growth and latex yield-related traits in Hevea brasiliensis Müll[END_REF][START_REF] Rosa | QTL detection for growth and latex production in a full-sib rubber tree population cultivated under suboptimal climate conditions[END_REF][START_REF] Souza | QTL mapping of growth-related traits in a full-sib family of rubber tree (Hevea brasiliensis) evaluated in a sub-tropical climate[END_REF]. However, for complex traits under the control of a large number of genes with small effects, such as yield, the efficiency of marker assisted selection approaches based on QTLs is limited, because it overestimates the effect of the strong QTLs while weak QTLs are not detected [START_REF] Muranty | Potential for marker-assisted selection for forest tree breeding: lessons from 20 years of MAS in crops[END_REF].

Currently, Hevea breeding involves within-family two-stage phenotypic selection (PS) followed by large-scale agronomic evaluation (Figure 1,left). Although a large number, i.e. several thousand, of full-sibs can be evaluated in the first stage (seedling evaluation trial, SET, with nonreplicated individuals), selection for rubber production at this stage is not very accurate [START_REF] Bombonato | Rubber tree ortet-ramet genetic correlation and early selection efficiency to reduce rubber tree breeding cycle[END_REF][START_REF] Gnagne | Méthodologie de sélection de l'hévéa : évaluation des seedlings[END_REF]. The second stage consists in small-scale clone trials (SSCT, with each genotype replicated in the form of several budded trees). The SSCTs make it possible to accurately assess clone yield, but the number of clones that can be evaluated in these trials is relatively low (< 200). This is followed by a long period of agronomic evaluation of growth rate, latex production, disease resistance, and other characteristics at the scale of tapped stands, in multi-local large-scale clone trials (LSCT).

GS is a very promising way to increase the rate of genetic progress in perennial crops [START_REF] Grattapaglia | Status and Perspectives of Genomic Selection in Forest Tree Breeding[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] because it allows the genetic value of a large number of selection candidates to be estimated at an early stage. In Hevea, the current SETs prior to clone trials could thus be replaced by more accurate genomic preselection. If GS is sufficiently accurate, it could even replace SSCTs. However, the decision to shift from a conventional PS scheme to a GS alternative calls for detailed studies. Indeed, the relative rate of genetic gain of different breeding approaches depends on their respective selection accuracy, selection intensity, and generation interval, with a trade-off between these parameters due to practical and economic constraints in the breeding program and to biological constraints of the species. Standard statistical methods for GS predictions include random regression best linear unbiased predictor (RR-BLUP) [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF], Bayesian least absolute shrinkage and selection operator regression (BLR) (de los [START_REF] De Los Campos | Predicting Quantitative Traits With Regression Models for Dense Molecular Markers and Pedigree[END_REF], and Bayesian reproducing kernel Hilbert Space (RKHS) [START_REF] Gianola | Reproducing Kernel Hilbert Spaces Regression Methods for Genomic Assisted Prediction of Quantitative Traits[END_REF].

BLR and RR-BLUP are linear approaches with different assumptions regarding the distribution of marker effects. Thus, RR-BLUP estimates marker effects following a normal distribution with common variance for all markers, while BLR uses a variance specific to each marker. RKHS is a semi-parametric and non-linear approach (i.e. using a non-linear genomic matrix) [START_REF] Pérez-Rodríguez | Comparison Between Linear and Non-parametric Regression Models for Genome-Enabled Prediction in Wheat[END_REF] that can capture both additive and non-additive effects [START_REF] Zhang | Establishment and Optimization of Genomic Selection to Accelerate the Domestication and Improvement of Intermediate Wheatgrass[END_REF]. When the purpose is to predict genetic values potentially including non-additive genetic effects (like clone values), it is appropriate to use models that take non-additive effects into account, either by modelling them explicitly (as it can be done in RR-BLUP and BLR) or implicitly (RKHS).

Despite the great economic importance of Hevea, no study has yet been published on the efficiency of GS compared with conventional PS in this species. Here, an alternative within-family breeding scheme for Hevea rubber production was suggested, in which the current phenotypic preselection of individual seedlings (SET) prior to clone trials (SSCT) would be replaced by genomic preselection in the nursery (Figure 1,right). As the GS model needs to be trained using phenotypic data, this alternative scheme would involve two SSCTs. The first, comprising a random sample of candidate clones (i.e. with no prior selection from SET results), would be used both to evaluate these candidates and to train the GS model; the second would be used to finalise the selection among the clones preselected by the GS model. The efficiency of the GS scheme compared with conventional PS will result from the accuracy and selection intensity (i.e. the number of clones genotyped to undergo genomic preselection) of GS. Genomic selection accuracy is usually estimated by within-site crossvalidation. However, such estimates may be biased upwards [START_REF] Beaulieu | Genomic selection accuracies within and between environments and small breeding groups in white spruce[END_REF]Lorenz et al., 2011, p.94;[START_REF] Ly | Relatedness and Genotype × Environment Interaction Affect Prediction Accuracies in Genomic Selection: A Study in Cassava[END_REF], and GS accuracy is consequently better estimated by validation using independent sites. The aim of this study was to carry out the first evaluation of GS for Hevea, using genotypic and phenotypic data on one family at two sites. For this purpose, two within-family clonal selection strategies for rubber production were compared: a new breeding scheme combining genomic preselection and PS, and the current conventional PS scheme. More precisely, (1) within-and between-site GS accuracy were estimated for rubber production in unobserved clones, (2) the effect of three parameters on GS accuracy was evaluated: statistical method of genomic prediction, size of training population, and molecular data (density and filtering), and (3) the increase in performance of the selected clones and in response to selection that could be expected from combining GS and PS compared with conventional PS was estimated using simulations based on the empirical data and on the between-site estimate of GS accuracy. Data on 330 clones from the F1 cross between two widely cultivated rubber clones (PB 260 × RRIM 600) were used, with phenotypic data collected from two independent clone trials in Côte d'Ivoire (189 clones at Site 1 and 143 clones at Site 2) and genomic data on 332 simple sequence repeat (SSR) markers.

Materials and methods

General overview

The study was divided into two parts. The aim of the first part was to obtain empirical estimates of GS accuracy for rubber production of unobserved clones of a F1 cross, with two independent field trials used for within-and between-site validations. The second part of the study aimed to estimate the additional annual response to selection that could be expected from combining GS and PS rather than using conventional PS. This was done by simulations based on the empirical data and on the GS accuracy estimated in the first part of the study.

For the first part (empirical estimation of GS accuracy), data on 330 clones from the F1 cross PB 260 × RRIM 600 were used. The clones were evaluated in two independent SSCTs in Côte d'Ivoire, with 189 clones at Site 1 and 143 at Site 2. The trials were implemented using conventional experimental designs, which allowed reliable estimations of clone values (hereafter referred to as phenotypes). The clones were also genotyped with 332 simple sequence repeat (SSR) markers. The GS model, trained using the molecular data and phenotypes of one part of the clones, predicted the phenotype of the other clones, for which molecular data only were used as inputs to the model. This made it possible to measure the accuracy of GS predictions of the performance of clones yet-to-be observed. The GS validation analyses were performed for predictions within and between sites to assess the usefulness of within-site accuracies (i.e. obtained by cross validation) for decision-taking regarding the practical implementation of GS. In addition, five standard statistical methods for genomic prediction were compared in terms of GS accuracy. Different training size, marker density, and SSR sampling method (sampling random SSRs or SSRs with the highest observed heterozygosity, Ho) were also used to quantify the effect of these three parameters on GS accuracy.

In the second part of the study (comparison of GS and PS schemes using simulation), the current conventional phenotypic breeding scheme (Figure 1, left) and an alternative scheme combining genomic preselection and PS (Figure 1, right) were simulated, and the two approaches were compared in terms of performance of the selected clones and of annual response to selection. The simulation was calibrated with the empirical data and with the results obtained in the first part of the study (genetic variance, PS accuracy, GS accuracy, etc.).

Empirical estimation of GS accuracy

Plant material and phenotyping

The parents of the F1 cross, PB 260 and RRIM 600, are two well-known and genetically unrelated clones that were selected in Malaysia. RRIM 600 originated from a cross made in 1937 (TJIR 1 × PB 86) and is the most widely planted clone in the world due to its high latex yield generated soon after tapping initiation and good adaptation to a variety of environments. Its potential for rubber production is medium. PB 260, issued from the cross PB 5/51 × PB 49, was obtained in 1958. It is a vigorous and high-yielding clone, largely used as female in crossings because it has one of the highest female fertilities among the best rubber clones used as parents, thus allowing for much larger progenies than other female parents [START_REF] Baudouin | Recurrent selection of tropical tree crops[END_REF]. It was recommended for plantation in Asia in the 1980s and 1990s and is still the second most widely planted clone in rubber-producing countries.

The two study sites are located in the coastal area of south-western Côte d'Ivoire: Site 1 (latitude: 4°40'54" N, longitude: 7°06'05" W, on the SOGB [Société des caoutchoucs de Grand Bereby] estate, elevation: 33 m a.s.l., gravelly clayey loam, with 189 clones), and Site 2 (latitude: 5° 19' 47.79" N, longitude: 4° 36' 39.74" W, on the SAPH [Société Africaine de Plantations d 'Hévéas] estate, elevation: 89 m a.s.l. deep sandy soil, with 143 clones). The two sites lie approximately 300 km from each other. The sites have a similar tropical climate, with 1,600 mm mean annual rainfall and a mean annual temperature of 26°C. Two clones were used at both sites, giving a total of 330 clones.

The two shared clones were not used for GS validation, but only to train the GS model. There was no preselection of the clones before the trials (i.e. no SET evaluation).

Site 1 and Site 2 trials were planted in July 2012 and July 2013 respectively, following almost complete block designs, with six blocks and individual trees randomised within each block, and with a planting density of 1,600 trees per hectare, with a spacing of 2.5 × 2.5 metres. Ramets were produced in the nursery by grafting on rootstocks generated from seeds issued from natural pollination of clone GT 1, and transplanted to the trials. The mean number of ramets per clone was 11 (range: 7 -17) at Site 1 and 13 (5 -20) at Site 2. This led to 2,016 ramets at Site 1 and 1,869 at Site 2.

Rubber production was recorded for each ramet according to the following protocol. The tapping system was in half-spiral on the trunk at 1 m above ground level, tapping every two days excepted on Sundays. In each trial, the six blocks were tapped by three tappers, with two blocks assigned to each tapper during the three consecutive months of the experiment. Each tapper tapped one block per day. No ethephon stimulation was applied to the trees in order to assess the natural latex flow of every tree. The latex was collected in plastic boxes of 180 ml attached to the trunks with rubber bands, and covers were screwed over the boxes between two tappings for preserving the latex production from rain. Every day in each trial, the boxes full of latex were collected and replaced by empty ones. The coagulated latex from the full boxes was extracted and manually pressed to eliminate the liquid serum. Tapping started 32 months after planting in Site 1 (end of dry season) and 38 months after planting in Site 2 (end of rainy season). For each ramet, the amount of rubber aggregated from all collected boxes during the 3-month tapping period (with a dry rubber content of around 65%) was computed. These raw production data were analysed for each site with a linear mixed model and the BLUP methodology, using the ASReml-R version 3.0 package [START_REF] Butler | Mixed models for S language environments: ASReml-R reference manual (Version 3)[END_REF]. This gave the clone genetic values, adjusted for effects related to the experimental designs (blocks) and for variations in size among the trees at the time of tapping (i.e. variations in girth of the trunk measured at 1 m above the ground just before initiation of tapping). These adjusted clone values are hereafter referred to as phenotypes. The broad sense heritability of clone mean level (H²) was calculated at each site as per Eq. (1).

= ( + ) (1)
Where is the genetic variance of the clones, the residual variance and ℎ the harmonic mean number of ramets per clone in the trial [START_REF] Gonçalves | Genetic variability for girth growth and rubber yield in Hevea brasiliensis[END_REF], with and obtained from the linear mixed model.

Marker genotyping

Leaf samples were collected on the original mother-trees of the clones issued from the seeds of the cross. Genomic DNA extraction and SSR genotyping were carried out following the method described by Le [START_REF] Le Guen | Genetic structure of Amazonian populations of Hevea brasiliensis is shaped by hydrographical network and isolation by distance[END_REF]. Site 1 clones were genotyped with 332 SSRs [START_REF] Tran | Genetic determinism of sensitivity to Corynespora cassiicola exudates in rubber tree (Hevea brasiliensis)[END_REF], and Site 2 clones were genotyped with a subset of 296 SSRs [START_REF] Achour | Cartographie génétique de la famille F1 PB260 x RRIM600 (Hevea brasiliensis) et détection de QTLs associés à la croissance[END_REF]. Table 1 lists the characteristics of the SSR molecular data obtained at each site. Sporadic missing SSR data 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], with parameters niterations set to 25 and nsamples to 20.

Statistical methods for genomic predictions

Three GS statistical methods were used to predict the genetic values of the validation clones: RR-BLUP, BLR and RKHS. In addition, BLR and RR-BLUP were carried out with two types of model, i.e. purely additive models (BLR_A and RR-BLUP-A) and additive plus dominance models (BLR_AD and RR-BLUP_AD). We did not consider the explicit modelling of epistatic effects for the sake of simplicity and assuming they would be negligible over additive and dominance effects.

For RR-BLUP_A and BLR_A, the model given by Eq. ( 2) was used.

y = 1µ + Zama + e (2)
Where y is the (k × 1) vector of phenotypes of training clones, k the number of clones, µ the overall phenotypic mean, 1 a column vector of 1s, ma the (n × 1) vector of allele additive effects, with n the total number of alleles, Za the (k × n) incidence matrix with elements = 0, 1 or 2 indicating the number of alleles j for clone i, and e the vector of residual effects following N(0, σ²e), with σ²e the residual variance. For RR-BLUP-AD and BLR_AD, the previous model was extended as indicated by Eq. ( 3).

y = 1µ + Zama + Zdmd + e (3)
Where md is the (p × 1) vector of dominance effects of all possible pairs of alleles at each SSR, with p the total number over all SSRs of possible combinations between two alleles of the same SSR, and Zd the (k × p) incidence matrix with elements = 1 or 0, indicating whether clone i possesses allele combination (pair) j or not. The genomic estimated genetic value (GEGV) of the validation clone i was obtained by Eq. ( 4) in RR-BLUP_A and BLR_A,

= ∑ (4)
and by Eq. ( 5) in RR-BLUP_AD and BLR_AD,

= ∑ + ∑ (5)
with the estimated additive effect of allele j, and the estimated dominance effect of the j th pair of alleles. For RR-BLUP, the ! " and ! # vectors were the BLUP solutions; and for BLR they were the posterior mean values over the post burn-in iterations. For BLR, σ²e followed a scaled inverse chi-square prior distribution, and ma and md followed conditional Gaussian prior distributions N(0, τ²aj σ²e) for allele j and N(0, τ²dj σ²e) for allele pair j, respectively. The τ²aj parameters were thus specific to each allele j, and the τ²dj to each allele pair j; and they followed exponential priors with rate λ²a / 2 and λ²d / 2, respectively, with the regularisation parameters λ²a and λ²d following gamma priors.

For RKHS, the model presented in Eq. ( 6) was used.

y = 1µ + g + e (6)
Where g = Kα is the vector of random genetic values of clones, K the (c × c) kernel constructed from the SSR data of the c clones, c the total number of clones (i.e. training and validation clones) and α the (c × 1) vector of regression coefficients to be inferred, with prior distribution N(0, Kσ²α). K gave the covariance structure among clones and had elements given by Eq. ( 7).

$ = % & '( (7)
With ) the squared Euclidean distance between clones i and j computed from their SSR genotypes, and h a bandwidth parameter. A multi-kernel approach based on a set of values of h (0.1, 0.5, 2.5) was implemented, as explained in Pérez and de los Campos (2013). The vector of GEGV * was obtained as per Eq. ( 8).

* = +, (8)
The BGLR R package version 1.0.5 [START_REF] Pérez | Genome-wide regression and prediction with the BGLR statistical package[END_REF] was used for BLR and RKHS with 30,000 iterations, with the first 9,000 as burn-in and a thinning interval of 10. For RR-BLUP, ASReml-R version 3.0 package [START_REF] Butler | Mixed models for S language environments: ASReml-R reference manual (Version 3)[END_REF] was used.

Validation approaches

The analyses were performed for predictions within and between sites, leading to four different validation approaches (Site 1 cross validation, Site 1 towards Site 2, Site 2 cross validation and Site 2

towards Site 1).

The clones from each site were randomly allocated into k sets used as validation replicates, with k=7 for Site 1 and k=5 for Site 2. In this way, the number of clones per set was similar in the validation experiments: 27 for Site 1 and 28 or 29 for Site 2, depending on the set. The allocation of clones to the validation sets was the same for all validation scenarios. Within-site validations were conducted using k-fold cross-validation approaches, successively using one of the k sets as the validation set and the remaining k-1 sets (or only some of them when varying the training size, see below) to train the GS model. For between-site validations, the k sets from one site (or some of the k sets when the training size was varied) were used to train the GS model, and the k sets of the other site were used for validation. The GS predictive ability was obtained for each set as the Pearson correlation between the GEGV (ĝ) and the phenotype (y) of the clones composing the set. Finally, GS accuracy was the predictive ability divided by the square root of the broad sense heritability (Lorenz et al., 2011, p. 94).

Effect of training size and SSR density

To quantify the effect of training size on GS accuracy, the number of sets composing the training population varied from one to k-1 (within-site validations) or k (between-site validations). For a given number of ns sets used for training, the different possible combinations of ns sets among the available sets (i.e. k-1 for within-site validations or k for between-site validations) were used successively.

These combinations of sets became training replicates. The resulting total number of replicates (validation replicates × training replicates) per validation experiment and training size varied from 5 to 175 (see Supplementary Table S. 1 for details).

To investigate the effect of SSR density on GS accuracy, different numbers of SSRs were also used, considering six levels of number of SSRs, from 10 to all SRRs. For a given number of SSRs, eight replicates of random samples of SSRs were made.

Effect of SSR sampling method

To investigate whether sampling SSRs with high observed heterozygosity (Ho) would lead to higher GS accuracy than randomly selected SSRs, Ho was computed for each SSR as the mean percentage of heterozygous individuals and the validations described above were run considering 12 levels of number of SSRs, from 10 to all SRRs. For a given number of SSRs, eight replicates of random samples of SSRs were made. With SSR sampling selecting the highest Ho, four replicates were also made, as some SSRs had the same Ho (in which case, the SSRs were chosen randomly). Here, all the clones were used to train the GS model.

Analysis of results

To study the effect of the statistical method for genomic prediction on GS accuracy, analyses of variance (ANOVA) were performed separately for each validation approach on the accuracy obtained using all SSRs and all the clones for training, with statistical method and validation replicate as factors. The mean levels of factors in the ANOVAs were compared using Tukey's honest significant difference test. To assess the effect of the SSR sampling method, the Wald-type permutation test of the R package GFD [START_REF] Friedrich | GFD: An R Package for the Analysis of General Factorial Designs[END_REF] was used, given that the assumptions of normality and variance homogeneity were not met, with SSR sampling method and validation replicate as factors.

The tests were carried out separately for each SSR level of each validation approach. Prior to these analyses, GS accuracy underwent Fisher's Z transformation.

Comparison of combined GS/PS and PS breeding schemes

The application of the conventional PS scheme on cross PB 260 × RRIM 600 (see Figure 1, left) was simulated, as well as a GS scheme in which clones of the same cross evaluated in a first SSCT would be used to train a GS model to make a preselection among unobserved clones (seedlings) of the same cross prior to their final evaluation in a second SSCT (see Figure 1, right). The simulation was calibrated with the results of the linear mixed model initially implemented to obtain the phenotypes, and with the results of the between-site empirical GS validations.

The simulation procedure started with the joint simulation of the true genetic values (TGV) (g), the seedling phenotypes in SET (y'), the estimated genetic values in SSCT (EGV, i.e. phenotypes) (y), and the genomic estimated genetic values (GEGV) (ĝ) of n individuals as per Eq. ( 9).

n = max(3,000, 190+nGS) (9)
With the 3,000 seedlings evaluated in SET, 190 clones evaluated in the first SSCT (used for both phenotypic selection and training of the GS model), and nGS the number of additional selection candidates allowed by GS (i.e. candidates subjected to genomic preselection at the nursery stage, prior to the second SSCT), with nGS varying from 100 to 5,000. These values were simulated using the mvrnorm function in the MASS R-package [START_REF] Venables | Modern Applied Statistics with S, Fourth[END_REF]. This required the variancecovariance matrix between g, y', y and ĝ given in Eq. ( 10 (10) and the mean phenotypic value of the clones µ, which were obtained as follows. The correlation between y' and y (< 76,7 ) was 0.34, taken from [START_REF] Gnagne | Méthodologie de sélection de l'hévéa : évaluation des seedlings[END_REF], and gave the correlation between rubber production in SET and SSCT for an unselected population related to the cross used in this study. For each site, the clone phenotypes given by the initial mixed model analyses were used to compute the associated variance, 7 . This initial analysis also gave µ, and the accuracy of phenotypic selection (< 0,7 , corresponding to the square root of ) in a SSCT with no preselection (genomic or based on SET evaluation). The variance of the TGV of the clones ( 0 ) was obtained as per Eq. ( 11) (Clark et al., 2012, appendix 1).

0 = 7 < 0,7 (11)
The GEGVs obtained from the empirical between-site validations were used to compute the associated variance, 0 . The GS accuracy (< 0,0 ) was taken from the between-site validations, and the correlation between y and ĝ was obtained as per Eq. ( 12) (Lorenz et al., 2011, p. 94;Muranty et al., 2015, appendix).

< 7,0 = < 0,0 < 0,7 (12)
Similarly, the variance of the seedling phenotypes ( 76) was calculated as per Eq. ( 13), the correlation between y' and g per Eq. ( 14), and the correlation between y' and ĝ per Eq. ( 15).

76 = 7 /< 7 > ,7 (13)
< 0,76 = 0 / 76 (14)

< 76,0 = 0 / 76 (15)
The mean values over the two sites were computed for each of these parameters and were used to calibrate the simulation.

The n simulated individuals served as starting point for the simulation of the conventional PS scheme and the alternative scheme combining GS and PS. For PS, a random set of 3,000 individuals was sampled among the n simulated individuals. Among them, the 190 individuals with the highest performance in SET (i.e. highest y') were retained to make the first SSCT, and the nsel = 10 clones with the highest EGV were finally selected among them. For combined GS/PS, a random set of 190 clones were sampled among the n simulated individuals to make the first SSCT. Then, the 185 clones with highest GEGV were selected among the nGS simulated clones (i.e. among those that were not evaluated in the first SSCT), to make the second SSCT (only 185 instead of 190 in the first SSCT, since in practice some clones from the first SSCT would be repeated in the second). Finally, nsel = 10 clones with the highest EBV were selected among the clones evaluated in the two SSCTs. The performance of the clones selected in combined GS/PS and PS schemes was computed as the mean TGV of the nsel selected clones. The annual selection response of the PS and combined GS/PS schemes was computed as the difference between the mean TGV of the nsel selected clones and the mean TGV of the n initial clones, divided by the number of years required to complete the breeding cycle (25 years) (Figure 1). The selection intensity of PS and combined GS/PS was computed as the mean EGV of the nsel selected clones and the mean EGV of the n initial clones, divided by the standard deviation of the EGV. The simulation process was repeated 5,000 times.

All analyses and simulations were conducted using the R software, version 3.4.1 (R Core Team, 2017).

Results

Phenotypic evaluations

Mean cumulated rubber production per tree was 78.7 g in Site 1 (range 0.50 -318.0) and 244.6 g in Site 2 (range 0.25 -840.1). The broad sense heritability of clone mean level (H²) was 0.9 at each site.

Statistical methods for genomic predictions

The GS accuracies obtained for rubber production were not affected by the statistical method used for predictions. When training the GS models with all clones and using all the SSRs, the mean GS accuracy over validation replicates ranged from 0.33 to 0.60 (Figure 2). However, this variation was mostly due to the validation approach, with statistical methods having a negligible effect. The differences in accuracy between statistical methods were thus not significant, regardless of marker density and size of training set, with no interaction found between the GS prediction method and SSR number (Supplementary Fig. Indeed, its mean accuracy across the four validation approaches when using all SSRs and all clones for training (0.498) was slightly higher than that of the other prediction methods (whose mean accuracy ranged from 0.488 for BLR_AD to 0.495 for RR-BLUP_A). Furthermore, it was the method that came out with the best average rank of the four validation approaches. When all the SSRs and all the clones were used for training, BLR_A gave a mean GS accuracy of 0.594 in Site 1 cross validation, 0.509 in Site 1 towards Site 2, 0.340 in Site 2 cross validation and 0.550 in Site 2 towards Site 1 validation.

Training population size and molecular marker data

GS accuracy for rubber production was strongly affected by the number of clones used to train the GS prediction model (training size) and by the number of SSRs (Figure 3).

GS accuracy increased with the training size regardless of validation approaches and number of SSRs used. For instance, when all the SSRs were used, increasing the training size from minimum to maximum values (i.e. by an average of 447.7%, from 296% in Site 2 cross validation to 600% in Site 1 towards Site 2), GS accuracy approximately doubled (mean of +93.6% across validation approaches, from 72.8% in Site 1 cross validation to 111.7% in Site 1 towards Site 2). With all validation approaches and numbers of SSRs, the increase in GS accuracy associated with increased training size followed a diminishing returns pattern. Thus, when 296 SSRs were used, increasing the training size from 28 to 56 clones increased GS accuracy by an average of 36.9% in the four validation approaches, while doubling the training size again to reach 111 clones increased the GS accuracy by slightly less (32.1%). Although usually GS accuracy did not reach a plateau, the shape of the curves showed that further increases in training sizes would have led only to minor additional gains in GS accuracy (except for Site 2 cross validation, due to the smaller overall population size). Similarly, with all the validation approaches and training sizes, GS accuracy increased with the number of SSRs.

Thus, increasing the number of SSR from minimum to maximum values (i.e. by an average of 3.0%, with 3.2% in Site 1 cross validation and 2.9% in other validations) when using the maximum training sizes, the average GS accuracy over validation approaches increased by 201.6% (from 134% in Site 1 cross validation to 296.2% in Site 2 cross validation). Again, a diminishing returns trend was observed for all validation approaches. For instance, with the largest training sizes, using 50 SSRs instead of 25 SSRs increased GS accuracy by 36.1% on average across all the validation approaches, while doubling again SSR density increased GS accuracy by 16.5% only. In Site 1 cross validation, for which more SSRs were available, using 332 SSRs resulted in the same accuracy as using 296. This indicated that no extra gain could be expected here from using more SSRs.

When SSR density was reduced, using the SSRs with the highest observed heterozygosity (Ho) generally resulted in significantly higher GS accuracies than using random SSRs (Figure 4). In particular, when the 125 to 200 SSRs with the highest Ho were used, GS accuracies were always significantly higher than the accuracies obtained with random SSRs, with an average increase of 13.9% (from 4.6% in Site 2 cross validation with 200 SSRs, to 21.1% in Site 2 towards Site 1 validation with 150 SSRs). Furthermore, in this range of number of SSRs, the Ho sampling approach led to almost always higher accuracies than using all SSRs, with an increase in GS reaching an average of 4.3% for the four validation approaches compared with using all the SSRs (the only exceptions being with 150 SSRs in Site 1 cross validation and with 200 SSRs in Site 2 cross validation, when GS accuracy with Ho SSR sampling was very slightly lower than when all SSRs were used). Mean GS accuracy of between-site validations thus reached 0.561, versus 0.530 using all SSRs. As expected, due to high variations in Ho among SSRs (Table 1), the SSR samples based on this parameter had a much higher mean Ho than the whole set of markers (Supplementary Fig. S 3).

Thus, when using 125 to 200 SSRs, mean Ho was 0.78, as against 0.64 with all the SSRs.

Validation approach

The effect of the validation approach on GS accuracy was investigated by comparing accuracies among validation approaches using the same training size and number of SSRs. In this case, withinlocation analysis gave much higher accuracies for Site 1 than for Site 2 (Figure 3). For instance, using 296 SSRs, within-Site 1 GS accuracy was 0.54 with 108 clones for training, versus only 0.34 for within-Site 2 accuracy with 115 training clones. By contrast, between-locations accuracies were similar when making predictions from Site 1 towards Site 2 and from Site 2 towards Site 1; and between-location GS accuracies were intermediate between the two within-site accuracies.

Site 1 cross-validation accuracy overestimated Site 1 towards Site 2 accuracy for all training sizes and numbers of SSRs (Figure 3). Thus, when using all the clones for training and all the SSRs, Site 1 cross-validation accuracy was 0.60, while Site 1 towards Site 2 accuracy fell to 0.51 (-14.9%).

By contrast, Site 2 cross-validation accuracy largely underestimated Site 2 towards Site 1 accuracy, for all training sizes and numbers of SSRs. Thus, when using all the clones for training and all the SSRs, Site 2 cross-validation accuracy was 0.34, while Site 2 towards Site 1 accuracy reached 0.54 (+61.7%).

Regarding the advantage of the Ho SSR sampling method over random sampling (Figure 4), consistent results were obtained between each within-site experiment and the between-site experiment in which the considered site was used for training and the other site for validation: in both cases, SSR sampling based on Ho gave higher accuracies than random sampling. In addition, the number of SSRs that gave the highest accuracy with Ho SSR sampling was the same in Site 2 cross validation and in Site 2 towards Site 1 validation (150), in both cases leading to higher GS accuracy than when all the SSRs were used. Similarly, although the number of SSRs that produced the highest accuracy with Ho SSR sampling in Site 1 cross validation and in Site 1 towards Site 2 validation differed (200 and 125, respectively), using the number of SSRs that gave the highest accuracy in Site 1 cross validation for Site 1 towards Site 2 validation would still have increased GS accuracy compared with using all SSRs.

Comparison of combined GS/PS and PS breeding schemes

The variance-covariance matrix between g, y, and ĝ used to calibrate the simulation is given Figure 5.

The mean phenotypic value (aggregated amount of rubber) was 186 g. The GS accuracy (< 0,0 ) was 0.561, corresponding to the mean accuracy obtained in between-site validations with the 125 to 200

SSRs with the highest Ho. The accuracy of SET (< 0,76 ) was 0.358.

The simulation showed that combining GS and PS outperformed conventional PS in terms of rubber production of the selected clones and annual selection response when genomic preselection was applied to a sufficient number of candidates, i.e. at least 1,000. In this case, additional rubber production was observed in the clones selected using GS (Figure 6). With 1,000 candidates, this additional production was very low (+0.4%) but increased when more candidates were used for preselection, and reached 5.9% when preselection was applied to 5,000 candidates. This led to an increase in annual response to selection when combining GS and PS compared with conventional PS, which started from +1% with 1,000 clones subjected to genomic preselection and reached +15% with 5,000 candidates for genomic preselection (Figure 7). The results also indicated that using a larger population of candidates for genomic preselection would have further increased the superiority of combined GS/PS over conventional PS, albeit only slightly. Genotyping 3,000 candidates for genomic preselection appeared as a good compromise between genotyping effort and efficiency of combined GS/PS (+4.2% mean production for the selected clones, corresponding to a +10.3% increase in annual selection response). In contrast, combining GS and PS performed worse than conventional PS when 100 and 500 candidates only were used for the genomic preselection.

As these values were means of 5,000 replicates of the simulation, they show the average extra gain that would result from the application of genomic preselection in many replicates of the F1 cross studied here. This is of major interest for breeders, but the actual gain that would be achieved in a given breeding program is also crucial. To assess this (and in particular to assess the probability that a given application of the GS scheme indeed performs better than the current PS scheme), Figure 7 also shows the distribution of the relative performance of combining GS and PS compared with conventional PS, in the form of a boxplot for each number of candidates for genomic preselection. For instance, the figure shows that with 2,000 candidates for genomic preselection, although the mean expected extra annual response to selection generated by GS reaches 7%, the first quartile is only slightly above the value corresponding to a similar performance by PS. This indicates that, although on average over a large number of replicates combining GS and PS will be better than conventional PS, for a specific application there is an almost 25% risk that GS would actually not perform better, or even worse, than PS (with the lowest value obtained being an annual selection response of GS reaching only 78.3% that of PS). Therefore, the best way to decide on the size of the population of selection candidates for genomic preselection is to consider both the mean expected annual selection response of GS and the distribution of the possible values around this mean. Thus, using 3,000 candidates, 75% of the simulation replicates gave an annual selection response of combined GS/PS at least 4.5% higher than when using conventional PS (with the maximum value reaching +45.6%), and the risk of GS actually performing worse than PS was low, at 9.2%. When increasing the number of candidates to 5,000, this risk dropped to 4.6%.

The increase in the relative performance of combined GS/PS compared with conventional PS when more candidates are used for genomic preselection resulted from the associated increase in selection intensity of combined GS/PS. Selection intensity in combined GS/PS was 15% lower than PS when 100 candidates were used for genomic preselection, but became roughly equivalent to PS with 1,000 candidates. It further increased to reach +14% when 5,000 candidates were used for the genomic preselection (see Supplementary Fig. S 4 for details). The fact that the selection intensity of the combined GS/PS scheme with 1,000 candidates to genomic preselection was similar to the selection intensity of the PS scheme with 3,000 individuals in SET resulted from the existence of two stages of selection. Indeed, genomic preselection can retain elite clones for SSCT that could have been discarded from the SET results, since the accuracy of genomic predictions is higher than SET accuracy. As a consequence, the 10 best clones selected at the end of the SSCT tend to perform better if the SSCT is preceded by genomic preselection rather than SET, leading to a higher selection differential and thus higher selection intensity in the combined GS/PS scheme than in the conventional PS scheme.

Finally, the better performance of the combined GS/PS scheme compared with conventional PS was the consequence of the greater selection accuracy of genomic preselection compared with phenotypic preselection with SET (GS accuracy being 56.7% higher) and of the greater selection intensity achieved when the number of candidates to genomic preselection was sufficiently high (≥1,000).

Discussion

The results presented here showed that applying the suggested breeding scheme combining GS and PS can increase rubber production in the cross PB 260 × RRIM 600. However, the advantage of this new breeding scheme over conventional PS resulted from GS accuracy, genetic variance and selection accuracy in SET and in SSCT, which vary among single crosses and traits. In particular, even in the case of GS implemented within full-sib families like here and despite the existence of deterministic equations, it remains difficult to predict GS accuracy for a particular trait in a given family [START_REF] Schopp | Genomic Prediction Within and Across Biparental Families: Means and Variances of Prediction Accuracy and Usefulness of Deterministic Equations[END_REF]. The study therefore needs to be extended to other families and traits, in particular using contrasted F1 crosses in terms of genetic and phenotypic variation.

Relevance of within-family GS for Hevea

The within-family GS scheme investigated here will not require restructuring breeding activities, already organised around full-sib families, and this is clearly a practical advantage for breeders. In addition, breeding schemes in which selection is applied within single crosses (i.e. full-sib families) are favourable situations for GS. In such biparental populations, there is a high linkage disequilibrium between marker alleles and gene alleles, which reduces the required marker density (as full-sibs share large chromosome segments), and there is no group structure [START_REF] Crossa | Genomic Selection in Plant Breeding: Methods, Models, and Perspectives[END_REF][START_REF] Lin | Genomic selection in crops, trees and forages: a review[END_REF].

Good results of within-family GS as implemented here have been reported in other plant species, with GS accuracies reaching moderate (i.e. between 0.5 and 0.7, as in the present study) to high values. For instance, GS accuracy estimated with a single-site cross validation was around 0.6 in a family of 500

Sitka spruce clones [START_REF] Fuentes-Utrilla | QTL analysis and genomic selection using RADseq derived markers in Sitka spruce: the potential utility of within family data[END_REF] and between 0.59 and 0.91 in a family of 180 Citrus clones [START_REF] Gois | Genome wide selection in Citrus breeding[END_REF].

A possible drawback of the within-family GS approach presented here is that it might not always be possible to obtain a training population of sufficient size. Hevea breeding programs use several families with limited resources, and the size of each family is therefore constrained. With the family used here, it appeared that using 175 clones to train the GS model was enough. However, this figure is close to the maximum amount of resources breeders can invest in a single family, and some families could require a larger training size, depending on their level of genetic variation. An alternative to the within-family approach studied here could be to implement GS in a population comprising several interconnected families, obtained using incomplete diallel or factorial mating designs. Although such a population would not be easy to obtain in Hevea due to the species' low female fertility, a comparison with within-family GS would be informative. This type of GS approach is implemented in a number of perennial species, including loblolly pine, spruce, eucalyptus [START_REF] Grattapaglia | Status and Perspectives of Genomic Selection in Forest Tree Breeding[END_REF], apple (eg Kumar et al., 2015;[START_REF] Muranty | Accuracy and responses of genomic selection on key traits in apple breeding[END_REF], and citrus [START_REF] Minamikawa | Genome-wide association study and genomic prediction in citrus: Potential of genomics-assisted breeding for fruit quality traits[END_REF]. This is interesting as it leads to a single (and therefore larger) training population compared with the various family-specific training populations required for the within-family GS approach. However, this increase in training size, although beneficial for GS accuracy, would be offset by a decrease in relatedness between the training set and the application set, a situation known to have a negative impact on GS accuracy. Therefore, in practice, a GS approach using a complex population involving several families could be more complicated to manage, with GS accuracy varying among selection candidates depending on their actual relationship with the training individuals. This could also actually lead to lower GS accuracies than family-specific training populations [START_REF] Crossa | Genomic Selection in Plant Breeding: Methods, Models, and Perspectives[END_REF][START_REF] Lenz | Factors affecting the accuracy of genomic selection for growth and wood quality traits in an advanced-breeding population of black spruce (Picea mariana)[END_REF][START_REF] Schopp | Genomic Prediction Within and Across Biparental Families: Means and Variances of Prediction Accuracy and Usefulness of Deterministic Equations[END_REF][START_REF] Toro | Accuracy of genomic withinfamily selection in aquaculture breeding programmes[END_REF][START_REF] Würschum | Accuracy of within-and among-family genomic prediction in triticale[END_REF]. In addition, from a practical point of view, the time needed to achieve and release a commercial clone could be longer with a complex multiparental population than with separate F1 families. This has to be taken into consideration as it represents a risk for Hevea breeding, where cycles are long and the resources invested in breeding activities are very limited.

Comparison of combined GS/PS and PS breeding schemes

The most important point for breeders regarding GS is the annual selection response that could result from its use, compared with the annual selection response of PS [START_REF] Resende | Assessing the expected response to genomic selection of individuals and families in Eucalyptus breeding with an additive-dominant model[END_REF].

Although PS and GS selection accuracies play a crucial role in this comparison, other factors that affect annual genetic gain must also be taken into consideration, i.e. relative generation interval and selection intensity of PS and GS. A few studies have ventured beyond estimating empirical GS accuracies and have used these estimates to evaluate the possible gain in annual selection response that GS could elicit. In eucalyptus, GS annual selection response is expected to be 50% to 300% greater than that of current PS, depending on the reduction in the duration of the breeding cycle and on GS selection intensity [START_REF] Resende | Genomic selection for growth and wood quality in Eucalyptus: capturing the missing heritability and accelerating breeding for complex traits in forest trees[END_REF][START_REF] Resende | Assessing the expected response to genomic selection of individuals and families in Eucalyptus breeding with an additive-dominant model[END_REF]. In black spruce, annual selection response should be 200% higher with the GS approach than with conventional selection, thanks to the shorter GS breeding cycle [START_REF] Lenz | Factors affecting the accuracy of genomic selection for growth and wood quality traits in an advanced-breeding population of black spruce (Picea mariana)[END_REF]. In Citrus, annual selection response is expected to increase by 31% to 420%, depending on how much the breeding cycle is shortened and on the trait concerned [START_REF] Gois | Genome wide selection in Citrus breeding[END_REF].

In Hevea, like in other perennial crops, the full potential of GS will be achieved over consecutive breeding cycles. Given the data available for this first GS study in this species, it was only possible to consider a single breeding cycle, whose duration could not be reduced due to the need for a SSCT to train the GS model. This explains why the increase in annual selection response reported here may seem modest compared with that reported in studies on other perennial crops. However, beyond the first cycle, breeding cycles will become shorter: only one SSCT will be required, since the GS model will have been calibrated with data from the first cycle. In addition, the training population used in the second cycle will comprise the aggregated data of the two SSCTs of the first cycle, and in the following cycles the data of the new SSCTs will be added to the training population. This is known to enhance GS accuracy [START_REF] Auinger | Model training across multiple breeding cycles significantly improves genomic prediction accuracy in rye (Secale cereale L.)[END_REF][START_REF] Cros | Training genomic selection models across several breeding cycles increases genetic gain in oil palm in silico study[END_REF][START_REF] Denis | Efficiency of genomic selection with models including dominance effect in the context of Eucalyptus breeding[END_REF]. Further studies are needed to investigate the efficiency of GS over several cycles in Hevea.

Another possibility would be to consider a GS scheme with only one SSCT in which the genomic predictions would be used to select clones in the second nursery, before their final evaluation in LSCT, instead of using GS to make a preselection before SSCT. This would have the advantage of reducing the duration of the breeding cycle. However, a simulation similar to the one presented here showed that, within the range of the number of selection candidates that can reasonably be genotyped, this approach was not advantageous in terms of annual selection response because the steep decline in accuracy between the SSCT (0.95) and the genomic predictions (0.561 on average over the two between-sites validations, i.e. a 40.9% decrease) was not offset by the shorter generation interval and/or higher selection intensity made possible by GS (data not shown). Our study therefore focused on a GS scheme in which the use of GS methodology was limited to the replacement of the conventional seedling evaluation trials prior to clone trials, and it showed this was sufficient to enhance the efficiency of the breeding scheme. A similar result was obtained in an oil palm study [START_REF] Cros | Genomic preselection with genotyping-by-sequencing increases performance of commercial oil palm hybrid crosses[END_REF], which evaluated the usefulness of genomic preselection prior to field evaluation, i.e. without reducing the breeding cycle duration, like in the present study. It thus showed that genomic preselection would increase bunch production by 6.5% to >10% when 2,000 to 10,000 candidates are used for genomic preselection.

Here, we used a single PS breeding scheme in order to benchmark the breeding scheme combining GS and PS. However, several PS schemes are possible. For instance, [START_REF] Gireesh | Fast track evaluation and selection of Hevea brasiliensis clones from a clonal nursery[END_REF] suggested the use of clonal nursery trials to optimize phenotypic breeding. It would therefore be interesting to implement new simulation studies to consider a broader range of possible PS and GS schemes.

Within-site and between-site accuracies

The between-site GS accuracies obtained in this study and the resulting estimate of annual selection response are appropriate for the environment considered here. The higher GS accuracy obtained in Site 1 cross validation compared with Site 1 to Site 2 validation was expected from the literature, which indicated that within-site cross-validations can lead to upward biases in GS accuracy [START_REF] Beaulieu | Genomic selection accuracies within and between environments and small breeding groups in white spruce[END_REF]Lorenz et al., 2011, p.94;[START_REF] Ly | Relatedness and Genotype × Environment Interaction Affect Prediction Accuracies in Genomic Selection: A Study in Cassava[END_REF]. For instance, in another perennial crop, black spruce, [START_REF] Lenz | Factors affecting the accuracy of genomic selection for growth and wood quality traits in an advanced-breeding population of black spruce (Picea mariana)[END_REF] obtained GS accuracy from between-site validation lower than accuracies obtained within the same site. The relatively small difference between the accuracies of Site 1 cross-validation and of Site 1 to Site 2 validation, and the similar accuracies when making predictions from Site 1 towards Site 2 and from Site 2 towards Site 1, indicated that genotype × environment (G × E) interactions, that could have been generated by differences in locations and years, were weak -probably because the two environments were similar. However, significant G × E interactions can occur in Hevea (see for example [START_REF] Da | Genotype-environment interaction and the number of test sites for the genetic improvement of rubber trees (Hevea) in São Paulo State, Brazil[END_REF][START_REF] Gonçalves | Genetic variability for girth growth and rubber yield in Hevea brasiliensis[END_REF][START_REF] Tan | Genotype x environment interaction studies in rubber (Hevea) clones[END_REF], and in this case the between-site GS accuracy would certainly be lower. In this case, the solution would be to take the environment into account in the prediction model. For this purpose, rubber geneticists will benefit from the methodology developed in cereals and legumes, where G × E modelling in the context of GS has been extensively studied [START_REF] Crossa | Genomic Selection in Plant Breeding: Methods, Models, and Perspectives[END_REF]. Surprisingly, GS accuracy obtained in Site 2 cross-validation was lower than the GS accuracy found in Site 2 towards Site 1 validation. What determined this result at this site remains unclear.

The effect of number of markers and SSR sampling method (random or based on high Ho)

observed for a single-site cross validation was in good agreement with the results obtained when a GS model was calibrated at this site to predict the values of clones evaluated at the other site. This indicated that, in the environment considered here, a single-site cross validation experiment made it possible to identify the number of SSRs and the method for choosing the SSRs that would yield the best GS accuracy that can be expected from using this experiment to train a GS model for predicting the rubber production of clones at another site.

Molecular data

In this study SSR markers were used, whereas in the vast majority of GS experiments in animals and plants, genotyping is carried out with single nucleotide polymorphism (SNP) markers.

Simple sequence repeats were used here because this type of marker has already been shown to be efficient in GS validation studies in oil palm [START_REF] Cros | Genomic selection prediction accuracy in a perennial crop: case study of oil palm (Elaeis guineensis Jacq.)[END_REF][START_REF] Marchal | Multivariate genomic model improves analysis of oil palm (Elaeis guineensis Jacq.) progeny tests[END_REF], table grapes [START_REF] Viana | Genome selection in fruit breeding: application to table grapes[END_REF], and flax [START_REF] You | Accuracy of genomic selection in biparental populations of flax (Linum usitatissimum L.)[END_REF]; and also because the biparental nature of the plant material used here suggested that the marker density achievable with SSRs could be sufficient. In the present study, 300 SSRs, which is a rather low density compared with what is usually found in GS studies, turned out to be sufficient to achieve the maximum GS accuracy that could be reached here.

However, this result holds for the F1 cross and for the training population size considered here, and it is possible that, in a different situation (for example with a larger training population), the GS accuracy would benefit from the use of more markers. Also, with the dataset considered here, it was possible to further reduce marker density with a slight increase, or at least no loss, in GS accuracy by using a subset of the 125 to 200 SSRs with the highest Ho. With multi-allelic markers in a single cross between heterozygous parents, Ho actually indicates how informative the markers are. Thus, the SSRs with Ho=1, which was the case for 25 to 50 SSRs per validation (Supplementary Fig. S 3), were those for which the two parents RRIM 600 and PB 260 had no alleles in common. When the two parents were heterozygotes, this corresponded to a situation with a balanced representation of the four alleles in the cross (the frequency of each allele being around 25%). This suggests that the marker density required to reach maximum GS accuracy is likely to vary among F1 crosses, depending on parental relatedness and heterozygosity. Other parameters were used for SSR screening (polymorphism information content (Botstein et al., 1980, p. 320) and expected heterozygosity, He) but preliminary analyses indicated that filtering using Ho yielded better results (data not shown).

The practical implementation of GS will require a high throughput and a cost-effective genotyping method to make the screening of large populations of selection candidates feasible. Even a reduced panel of SSRs might not be competitive in terms of cost compared with genotyping approaches involving SNPs. In addition, if the method is implemented over several generations, it will probably be necessary to increase marker density in order to limit decline in accuracy (Grattapaglia, 2017, p. 216). To our knowledge, there is currently no SNP array available in Hevea, but genotyping by sequencing (GBS) [START_REF] Elshire | A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species[END_REF], which has already been used in this species to construct a high density linkage map [START_REF] Pootakham | Construction of a high-density integrated genetic linkage map of rubber tree (Hevea brasiliensis) using genotyping-by-sequencing (GBS)[END_REF], could generate the molecular data required for GS in Hevea. Furthermore, approaches specific to biparental crosses that combine GBS and a relevant imputation methodology could be used to further increase the cost efficiency of large-scale genotyping [START_REF] Gorjanc | Potential of Low-Coverage Genotyping-by-Sequencing and Imputation for Cost-Effective Genomic Selection in Biparental Segregating Populations[END_REF][START_REF] Technow | Parent-progeny imputation from pooled samples for costefficient genotyping in plant breeding[END_REF].

Models and statistical methods for genomic predictions

Like in the present study, empirical GS accuracies have frequently been found to be unaffected by the statistical method of prediction (here BLR, BRR, and RKHS). Several examples with similar results are thus available in perennial crops. For various growth traits in eucalyptus, similar accuracies were obtained using BLR, RR-BLUP, and RKHS by Tan et al. (2017a), and using BLR and Bayesian random regression (BRR, similar to RR-BLUP) by [START_REF] Müller | Genomic prediction in contrast to a genome-wide association study in explaining heritable variation of complex growth traits in breeding populations of Eucalyptus[END_REF]. In oil palm, BLR and BRR gave the same accuracies for yield components [START_REF] Cros | Genomic selection prediction accuracy in a perennial crop: case study of oil palm (Elaeis guineensis Jacq.)[END_REF].

The explicit modelling of dominance effects with BLR_AD and RR-BLUP_AD had no effect on GS accuracy. Simulations in eucalyptus showed that including dominance in the GS model for prediction of clone performance improved accuracy when dominance effects were preponderant (ratio of dominance to additive variance of 1.0) and heritability was high (H²=0.600) [START_REF] Denis | Efficiency of genomic selection with models including dominance effect in the context of Eucalyptus breeding[END_REF]. Simulations in loblolly pine showed that including dominance in the GS prediction model improved accuracy when the ratio of dominance to phenotypic variance was over 20% [START_REF] De Almeida Filho | The contribution of dominance to phenotype prediction in a pine breeding and simulated population[END_REF]. With empirical data on eucalyptus, [START_REF] Tan | Genomic relationships reveal significant dominance effects for growth in hybrid Eucalyptus[END_REF] reported that GS accuracy for traits with large dominance variance was increased by including dominance effects in the model.

However, in apple, [START_REF] Kumar | Genome-Enabled Estimates of Additive and Non-additive Genetic Variances and Prediction of Apple Phenotypes Across Environments[END_REF], empirically obtained similar GS accuracies with models with or without non-additive effects for fruit quality traits, despite a high proportion of non-additive variance in some traits. This apparent discrepancy could come from the fact that [START_REF] Kumar | Genome-Enabled Estimates of Additive and Non-additive Genetic Variances and Prediction of Apple Phenotypes Across Environments[END_REF] used a training population of around 230 individuals, much smaller than that used by the previously cited authors, who used training populations of at least 800 individuals. It can therefore be hypothesised that, in the present study, including dominance effects in the GS models did not affect accuracy because dominance variance was not large enough and/or because the training populations were too small (from 114 to 189 individuals). Similar reasons are likely to explain the fact that RKHS did not perform better than the other methods.

Conclusions

The within-family GS strategy investigated here will lead to the release of more productive Hevea clones than clones selected with the current PS scheme. This will increase the yield of rubber from existing plantations, and thus help to meet the demand for natural rubber while minimising environmental costs. With a F1 cross between two widely cultivated clones, PB 260 × RRIM 600, a mean empirical GS accuracy of 0.53 was obtained in predictions between two independent sites when using all the clones for training and all the SSRs. SSR density and training size markedly affected GS accuracy. Mean between-site GS accuracy reached 0.561 when using the 125 to 200 SSRs with the highest Ho. In contrast, the statistical method used to obtain the genomic predictions of clone values did not affect GS accuracy. Based on this empirical result, simulations showed that by applying a genomic preselection among 3,000 seedlings in the nursery prior to clone trial, instead of the current low-accuracy phenotypic preselection on 3,000 seedlings, the rubber yield of the clones selected in the F1 cross considered would have been 4.2% higher, corresponding to a 10.3% increase in annual selection response. This resulted from the greater selection accuracy of genomic preselection compared with phenotypic preselection.

The results presented here showed that combining GS and PS can increase rubber production in the cross PB 260 × RRIM 600. However, before generalising GS in rubber breeding, this study needs to be extended to other families because the results obtained, and in particular the GS accuracies and selection response, are affected by the genetic characteristics of the parents of the F1 cross used.

Similarly, studies considering other traits, such as growth and architecture, are needed. It is also necessary to compare GS and PS in terms of selection response per unit cost and to investigate the efficiency of GS over consecutive breeding cycles, which will make it possible to shorten the breeding cycle in the cycles following model training. Furthermore, using a broader range of environments for between-site validations will be of major interest.

Tables Table 1 Characteristics

  S 1), nor between the GS prediction method and training size (Supplementary Fig.S 2). For the rest of the study, only the BLR_A GS prediction method was used.

FiguresFigure 1 .

 1 Figures

Figure 2 . 37 Figure 3 .

 2373 Figure 2. GS accuracy for rubber production according to statistical method of GS prediction, and

Figure 4 .

 4 Figure 4. GS accuracy in predicting rubber yield according to SNP sampling method (highest observed

Figure 7 .

 7 Figure 7. Annual response to selection in the GS scheme, expressed in % of the annual selection

  of the simple sequence repeat (SSR) molecular data obtained at each site. Ho:

	observed heterozygosity
		Site 1	Site 2
	Number of SSRs	332	296
	Missing data (%)	2.7%	2.2%
	Range of missing data (%) per SSR	0.0% -58.1%	0.0% -51.0%
	SSRs with ≤5% missing data (%)	87.3%	92.2%
	Range of missing data (%) per clone	0.0% -21.0%	0.0% -32.8%
	Clone with ≤5% missing data (%)	88.0%	95.1%
	Mean number of alleles per SSR (range)	2.56 (2 -4)	2.56 (2 -4)
	Total number of alleles	850	759
	Mean allele frequency (range)	0.39 (0.14 -0.86) 0.39 (0.15 -0.84)
	Mean Ho per SSR (range)	0.64 (0.34 -1)	0.64 (0.33 -1)

Acknowledgments

We thank the Institut Français du Caoutchouc (IFC), and the companies Michelin, SIPH, and SOCFIN for funding the part of this research devoted to Hevea rubber in the framework of the IFC project "Hevea varietal creation in West-Africa". We also thank the companies SOGB (Société Grand Bereby) and SAPH (Société Africaine de Plantations d'Hévéa) for their logistical assistance in the field experiments in Côte d'Ivoire. We acknowledge Laurence Dedieu for reviewing the paper, and the CETIC (African Center of Excellence in Information and Communication Technologies) for support.

Data availability

The datasets generated and analysed during the current study are available from the corresponding author.

Conflict of interests

The authors declare no conflict of interest.

Author contributions

DC carried out data analysis and wrote the paper, with help of ACD and VLG. LM, JO, and JB carried out preliminary data analysis. AM and MS provided assistance and logistics for trial setting up and phenotyping. VLG and ACD designed field experiments, supervised collection of phenotypic data and generation of molecular data by DMT and ZA.