S Ben-Sadoun

Renaud Rincent

J Auzanneau

François-Xavier X Oury

B Rolland

E Heumez

Catherine Ravel

Gilles Charmet

Sophie Bouchet
email: sophie.bouchet@inra.fr

Economical optimization of a breeding scheme by selective phenotyping of the calibration set in a multi

Keywords: Multi-trait selection, Trait-assisted selection, Selective phenotyping, Resource allocation optimization, Wheat baking quality, Coefficient of determination (CD)

Introduction

Thanks to exponential genotyping cost decrease, computing and data storage capacities improvement, methodological and statistical developments, genomic selection is becoming an essential tool to optimize breeding efficiency by predicting the genetic value of candidates for complex traits using molecular marker information (Whittaker et al. 2000;[START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF]. The use of genomic predictions at key steps of a breeding program can improve genetic gain per unit of time and cost, by skipping some field evaluation steps or optimizing the experimental design, i.e. minimizing the number of lines or replicates evaluated in one environment [START_REF] Heslot | Optimization of selective phenotyping and population design for genomic prediction[END_REF].

To optimize phenotyping cost, we can also take advantage of correlated (secondary) traits that are easier or cheaper to measure, or can be obtained at an earlier stage. It is possible to predict correlated traits simultaneously using multivariate best linear unbiased prediction (BLUP; [START_REF] Henderson | Multiple Trait Evaluation Using Relatives' Records[END_REF]. Multi-trait genomic prediction (MT) models benefit from information contained in both genetic correlation between traits and genetic relationship among individuals [START_REF] Calus | Accuracy of multi-trait genomic selection using different methods[END_REF]. Two strategies are possible: 1) either the training population is genotyped and phenotyped for both traits and the candidate population is genotyped but not phenotyped for any of the traits, or 2) the candidates can be partly phenotyped for the trait of interest or the secondary trait. The latter strategy is called trait-assisted genomic selection (TA) [START_REF] Fernandes | Efficiency of multi-trait, indirect, and trait-assisted genomic selection for improvement of biomass sorghum[END_REF]. Benefits of MT models over singletrait genomic prediction (ST) models were reported using simulated and empirical data [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] Hayashi | A Bayesian method and its variational approximation for prediction of genomic breeding values in multiple traits[END_REF][START_REF] Guo | Comparison of single-trait and multiple-trait genomic prediction models[END_REF]. As for single trait prediction, under a major QTL genetic architecture, [START_REF] Jia | Multiple-Trait Genomic Selection Methods Increase Genetic Value Prediction Accuracy[END_REF] found that Bayesian multivariate models (BayesA or BayesCπ) performed better than multi-trait GBLUP model. But for traits with polygenic genetic architecture, multi-trait GBLUP model was equal to the Bayesian multivariate models. MT models can however suffer from a high computational demand, time and some convergence problems [START_REF] Michel | Improving the baking quality of bread wheat by genomic selection in early generations[END_REF]. Obviously, genetic correlation between traits is a key factor determining the MT advantage over ST methods [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] Hayashi | A Bayesian method and its variational approximation for prediction of genomic breeding values in multiple traits[END_REF][START_REF] Guo | Comparison of single-trait and multiple-trait genomic prediction models[END_REF]. Although MT models improve the predictive ability when the targeted trait has a low heritability and the secondary trait has higher heritability, the advantage of MT models to predict high heritability traits is low [START_REF] Jia | Multiple-Trait Genomic Selection Methods Increase Genetic Value Prediction Accuracy[END_REF][START_REF] Hayashi | A Bayesian method and its variational approximation for prediction of genomic breeding values in multiple traits[END_REF][START_REF] Guo | Comparison of single-trait and multiple-trait genomic prediction models[END_REF].

Other factors also affect the prediction accuracy, including the size and the composition of the training set.

Assuming that the number of markers is sufficient, the accuracy of the calibration model strongly depends on congruency between the allelic composition represented in the training population to estimate marker effects, and the allelic composition of the candidates whose performance is to be predicted [START_REF] Habier | The impact of Genetic Relationship Information on Genome-Assisted Breeding Values[END_REF]). When the prediction uses unrelated populations to train the prediction equations, prediction accuracy actually becomes negligible [START_REF] Crossa | Genomic prediction in CIMMYT maize and wheat breeding programs[END_REF]. Methods based on minimizing the mean of the prediction error variance (PEVmean; [START_REF] Rincent | Maximizing the Reliability of Genomic Selection by Optimizing the Calibration Set of Reference Individuals: Comparison of Methods in Two Diverse Groups of Maize Inbreds (Zea mays L.)[END_REF][START_REF] Akdemir | Optimization of genomic selection training populations with a genetic algorithm[END_REF][START_REF] Isidro | Training set optimization under population structure in genomic selection[END_REF][START_REF] Sarinelli | Training population selection and use of fixed effects to optimize genomic predictions in a historical USA winter wheat panel[END_REF] or maximizing the mean of the generalized coefficient of determination (CDmean) of the estimated contrast between candidates and the mean of the panel were proposed to optimize the sampling of key individuals to be phenotyped in unstructured [START_REF] Rincent | Maximizing the Reliability of Genomic Selection by Optimizing the Calibration Set of Reference Individuals: Comparison of Methods in Two Diverse Groups of Maize Inbreds (Zea mays L.)[END_REF]) and in structured multi-familial populations [START_REF] Rincent | Predicting genomic selection efficiency to optimize calibration set and to assess prediction accuracy in highly structured populations[END_REF]. In any case, accurate predictions will not be obtained until a large enough pool of individuals that is relevant to the candidates is genotyped and phenotyped with a sufficient precision.

In bread wheat breeding programs, baking tests or dough rheology measurements are part of official variety evaluation for registration in many countries. Improving baking quality is a key aspect of wheat improvement but criteria may be specific to countries and targeted markets. In France, the most important criterion for baking quality is a normalized (NF V03-716) bread-making score (BMS). It is used to classify registered varieties into end-use quality classes. Varieties presenting high BMS scores are reserved for French baking; they must yield at least 2% higher than reference varieties in registration trials. Varieties with intermediate BMS scores are mainly meant to be exported or used in mixed flours and must yield 4% higher than reference varieties. Varieties with lower BMS can only be used as feed or fuel and must yield 7% more than reference varieties. This BMS score is a complex quantitative trait that integrates several physical measurements of dough and loaf. Its assessment requires a large amount of grains and flour which are not available in the first generations of the breeding program, i.e. in rownursery or small plots. In addition, BMS evaluation is time-consuming and quite expensive. Therefore, this trait is only evaluated at the very last steps of wheat breeding programs after the number of candidate lines was largely reduced. In French breeding programs, BMS is often estimated at earlier stages using less expensive traits such as the Chopin alveograph parameters which describe the dough rheological characteristics. Dough strength (W), in particular has a correlation of approximately 0.4 to BMS (Oury et al. 2010).

Furthermore, dough strength was shown to be itself influenced by high-molecular-weight glutenin subunits (HMW-GS; reviewed by Shewry 2009) that are cheaper to type and necessitate just a few grains. HMW-GS are encoded by the Glu-1 loci, named Glu-A1, Glu-B1 and Glu-D1, located on the long arm of the homoeologous chromosomes of group 1 [START_REF] Payne | Genetics of Wheat Storage Proteins and the Effect of Allelic Variation on Bread-Making Quality[END_REF]. Each locus comprises the two tightly linked genes Glu-1-1 and Glu-1-2.

To predict the degree of elasticity or tenacity of the dough of candidates to selection, HMW-GS used to be typed by sodium-dodecyl-sulphate polyacrylamide gel electrophoresis (SDS-PAGE) of grain proteins. Several attempts were made to derive molecular markers from glutenin DNA sequences (reviewed by [START_REF] Liu | Functional markers in wheat: current status and future prospects[END_REF]) because they are cheaper and easier to implement routinely. [START_REF] Ravel | SNP markers for early identification of high molecular weight glutenin subunits (HMW-GSs) in bread wheat[END_REF] recently developed KASP markers that predict the rheological properties of the dough slightly better than SDS-PAGE.

Previous studies concerning multi-trait genomic predictions of baking quality related parameters such as dough rheological traits [START_REF] Hayes | Accelerating wheat breeding for end-use quality with multi-trait genomic predictions incorporating near infrared and nuclear magnetic resonance-derived phenotypes[END_REF][START_REF] Michel | Improving the baking quality of bread wheat by genomic selection in early generations[END_REF][START_REF] Lado | Resource allocation optimization with multi-trait genomic prediction for bread wheat (Triticum aestivum L.) baking quality[END_REF] did not focus on BMS. In addition, none of these approaches tested the economic gain achievable with MT and TA models.

In this study, we demonstrate how to decrease the cost of breeding an expensive trait by using the information of a cheap correlated trait. To do so, we used genomic prediction of BMS as a case study. We evaluated different strategies in terms of predictive ability and cost. We compared ST models, univariate genomic prediction models with molecular markers associated with HMW-GSs as fixed effects (ST-glu), MT and TA models using dough strength (W) as a correlated trait. Only individuals of the training set were phenotyped for the trait of interest in both MT and TA models. In MT model, only the individuals of the training set were phenotyped for the correlated trait W. In TA model, the training set and a variable proportion of the validation set were phenotyped for the correlated trait. In addition, we adapted the CDmean criterion [START_REF] Rincent | Maximizing the Reliability of Genomic Selection by Optimizing the Calibration Set of Reference Individuals: Comparison of Methods in Two Diverse Groups of Maize Inbreds (Zea mays L.)[END_REF]) to a multi-trait context (CDmulti) in order to determine which individuals should be phenotyped for the secondary and the trait of interest to optimize the predictive ability of the trait of interest.

Materials and Methods

Plant material

We analyzed a breeding population of F8-F9 winter-type bread wheat lines developed by the Institut National de la Recherche en Agriculture, Alimentation et Environnement (INRAE) and its subsidiary company Agri-Obtentions. Each of the 1912 lines were evaluated between 2000 and 2016. Each year, the lines were phenotyped in 7 to 9 locations in France (Supplementary Table S1). In each environment, some registered varieties were used as controls, (6.8 control lines on average for trials in which BMS was phenotyped). Each control line has been evaluated for BMS for several consecutive years, up to 8 years for some of them. Those controls are essential to estimate adjusted means in further analyses. Crop management methods corresponded to high yield objectives (optimized pesticide, fungicide and nitrogen amount) to avoid confusing effects of disease or N-supply limitation.

Genotyping data

In total, 814 lines were genotyped using a 280K SNP array [START_REF] Rimbert | High throughput SNP discovery and genotyping in hexaploid wheat[END_REF]. As predictive ability was the same using a subset of 35K markers that minimize Linkage Desequilibrum (LD; result not shown), the 57 breeding lines that were produced after 2015 were genotyped with a 35K SNP array included in the 280K array. Finally, 21210 markers with minor allele frequency superior to 0.05 and less than 10% of missing data were kept for analyses. Missing data were imputed using a kinship-based EM algorithm [START_REF] Poland | Genomic selection in wheat breeding using genotyping-by-sequencing[END_REF], based on the assumption that marker genotypes follow a multivariate normal (MVN) distribution. A genomic relationship matrix K was computed using Endelman and Jannink equation (2012):

K= 𝑊𝑊 𝑇 2 ∑(𝑝 𝑘 -1)𝑝 𝑘 (1)
where W is a centered N × M marker matrix whose i-th row and k-th column is wik = xik + 1 -2pk with Xik the genotype of the i-th individual for the k-th marker as {-1,0,1} and pk the allele frequency at the k-th marker. These algorithms were implemented in the A.mat function of the rrBLUP package [START_REF] Endelman | Ridge regression and other kernels for genomic selection with R package rrBLUP[END_REF].

Additionally, a subset of 200 lines was genotyped at 11 DNA-based KASP markers derived from the HMW-GS loci Glu-A1, Glu-B1, and Glu-D1 [START_REF] Ravel | SNP markers for early identification of high molecular weight glutenin subunits (HMW-GSs) in bread wheat[END_REF]. We calculated Pearson's correlation between each of the KASP marker and the 21210 SNP markers. We identified 18 SNP markers with correlation higher than 0.9 with at least one of the HMW-GS markers. In the model using those KASP markers as fixed effects, these 18 SNP markers were filtered out to calculate the kinship matrix.

The datasets with the genotyping data are available in the INRAE Dataverse repository (https://data.inra.fr/). They can be accessed with the following link https://doi.org/10.15454/AABGO7.

Phenotypic data

In total, 871 lines were both genotyped and evaluated for grain yield in 9.9 environments on average. Although yield, height, heading date and protein were recorded for all the F8 lines, baking and rheological traits were only available in half of the lines that were selected for F9 evaluation. Protein content was estimated by Near-Infrared-Reflectance. This test was realized according to American Association of Cereal Chemists (AACC 39-10) method, using a wholemeal flour produced on a Cyclotec mill with a 0.8 mm sieve. Chopin alveograph (method AACC 54-30A) was carried out on flour obtained with a Chopin mill or a Brabender Senior mill (the extraction rate of these two mills is about 70%). The dough strength (W), the tenacity (P) and the extensibility (L) of the dough were evaluated with a Chopin Alveolink apparatus (Tripette & Renaud, 92396 Villeneuve-la-Garenne, France). The bread-baking tests were realized according to the AFNOR method NF V03-716 (also called BIPEA test). This baking test gives a measure of loaf volume, a score for dough (out of 100 units), a score for crust (out of 100 units) and a score for crumb (out of 100 units); the sum of these 3 scores giving the final BIPEA score (out of 300 units).

This final BIPEA score corresponds to the Bread Making Score (BMS). In total, 398 lines were both genotyped using an array and phenotyped for BIPEA BMS.

The dataset with the phenotypic data is available in the INRAE Dataverse repository (https://data.inra.fr/). It can be accessed with the following link https://doi.org/10.15454/AABGO7.

Statistical analysis of phenotypic data

Adjusted means of the lines were estimated using Best Linear Unbiased Estimations (BLUEs) in ASReml-R library [START_REF] Butler | ASReml-R reference manual[END_REF]) using the following model for each trait:

yijk = μ + gi + ej + gi × ej + εijk (2)
where yijk is the phenotypic value of the i-th genotype in the j-th environment (combination of year and site) and for the k-th repetition, μ is the overall mean, gi is the fixed effect of the i-th genotype, e j is the effect of j-th environment, gi × ej is the effect of the interaction between the i-th genotype and the j-th environment, and εijk is the residual error for the i-th genotype in the j-th environment and for the k-th repetition. All effects except the genotype effect were considered as random. Correlation between traits was calculated as the Pearson's correlation between BLUEs. Note however that traits linked to baking quality were measured only once in each environment.

Therefore, the gi × ej interaction was not included in the model for these traits.

The estimation of adjusted means can affect the predictive ability of the genomic prediction models. Therefore, we evaluated the impact of the number of phenotypic records used to calculate BLUEs on the predictive ability.

We especially investigated the importance of control lines that were phenotyped in all locations during several successive years in the estimation of BLUEs and the predictive ability, by comparing two different cases of BLUE calculations:

1) In the first case, all phenotypic records for control lines and one to four phenotypic records (in four different trials) for the other lines were used to estimate BLUEs for BMS and W;

2) In the second case, one to four phenotypic records for each line (the control and the other lines) were used to estimate BLUEs for BMS and W.

For each scenario, the phenotypic records were randomly sampled 10 times.

The variance components of the traits were assessed with ASReml-R using the model (2) but in this case all the effects (including genotypic effects) were considered as random.

Repeatability was estimated at the plot and design levels as:

plot repeatability = 𝜎 𝐺 2 𝜎 𝐺 2 +𝜎 𝐺×𝐸 2 +𝜎 ε 2 (3) design repeatability = 𝜎 𝐺 2 𝜎 𝐺 2 +𝜎 𝐺×𝐸 2 𝑡 ⁄ +𝜎 𝜀 2 (𝑡×𝑟) ⁄ (4)
Where σ 2 G; σ 2 G×E and σ 2 ε are the genotype, genotype × environment and residual variances respectively, t and r are the average number of trials per genotype and replicates in each trial per genotype, respectively. The variance components and repeatability were estimated using lines which were both phenotyped and genotyped. For the traits linked to rheological quality, σ 2 G×E term was removed from equations ( 3) and (4).

Genomic prediction models

Single-trait GBLUP models

First, a single-trait GBLUP model was used to obtain genomic predictions for different traits:

y = Xβ + Za + ε (5)
Where y is a N × 1 vector of BLUEs obtained in the first step with N the number of lines; β is the vector of fixed effects; X is the corresponding design matrix for the fixed effect; Z is the design matrix for the random effects; a is the N x 1 vector of lines random effects with the additive genetic variance σ 2 a and a~N(0,Kσ 2 a); and ε is the residual error with the residual variance σ 2 ε and ε~N(0,Iσ 2 ε).

Three types of single-trait models were performed. The first one was a simple ST model in which the vector of fixed effects β modelled the grand mean. The second model, called ST-glu, included the 11 HMW-GS markers as fixed effect. The third model (ST-glu_stepwise) included as fixed effect the HMW-GS markers that explained a significant genetic variance. They were selected using a stepwise regression linear model implemented in the stepAIC function from the MASS package [START_REF] Venables | Modern Applied Statistics with S. Fourth Edition[END_REF]. In this model, no line random effect was included.

Multi-trait GBLUP models

The same genomic relationship matrix K was used for fitting the multi-trait model:

y = Xβ + Za + ε (6)
Where y corresponds to an N × 2 vector of BLUEs for two traits, BMS which is the trait of interest and W which is a correlated trait less expensive to measure, β is the vector of fixed effects (the grand mean in this study); X is the corresponding design matrix for the fixed effect; a is the vector of N × 2 line effects with the corresponding random effect design matrix Z and a∼MVN (0, Σa⊗K) with the variance-covariance matrix Σa of the form:

( 𝜎 𝑎1 2 𝜎 𝑎12 𝜎 𝑎12 𝜎 𝑎2 2 ) (7)
where σa1 2 and σa2 2 are the genetic variance of the first and second trait, respectively, and σa12 is the genetic covariance between both traits. The variance of the residual effect followed ε∼MVN (0, Σε⊗I) where I is the identity matrix of dimension N × N and Σε has the same form than Σa. ⊗ indicate the Kronecker product operator between matrices.

We investigated two ways of predicting the value of the trait of interest using multivariate models: either we predicted it for un-phenotyped individuals (MT model), or we predicted it for a set of individuals that were fully or partly phenotyped for the correlated trait (TA model). In both MT and TA scenarios, measurements of the targeted trait were only available for the training population.

Model validation

To compare the different models, predictions were evaluated using 5-fold cross-validations. This method consists in randomly splitting the data set in 5 folds of equal size and using 4 folds as the training set in order to predict performance of the lines from the fifth fold. Predictive ability is calculated for each fold as the Pearson's correlation between predicted values and BLUEs. The same procedure is repeated 20 times. Predictive abilities and their standard deviations were calculated as the mean of the 100 Pearson's correlations. Models were tested using the lines which were both genotyped and phenotyped and which were not considered as control lines.

Cost evaluation

Budget allocated to phenotyping was obtained using the following equation:

B = CpBMS × NBMS + CpW × NW (8)
Where NBMS and NW are the number of lines phenotyped for BMS and for W respectively; and CpBMS (150€) and Cpw (20€) are the cost of phenotyping one line for BMS and for W respectively.

In the equation ( 8) CpBMS and Cpw are numeric constants but NBMS and NW are variable depending on the experimental design. It means that for a fixed budget there is several ways to allocate it between lines that are phenotyped for BMS or W, in the training set and in the validation set.

Phenotyping cost for W corresponds to a minimum cost when we do it internally at INRA. BMS must be evaluated by certified companies and the price can vary between providers and depending on the volume we phenotype each year. In order to evaluate the impact of the cost ratio between W and BMS on the predictive ability of BMS, we made this ratio vary (ratio = 1/2, 1/3, 1/4, 1/5, 1/10 and 1/20). We performed this analysis with 3 different budgets, with CpBMS=150€, with a MT approach and a TA approach where all lines are phenotyped for W (Supplementary Table S2).

Multi-trait CDmean (CDmulti) and optimization algorithm

The objective was to optimize the choice of individuals from the training set that should be phenotyped for the trait of interest BMS and for the correlated trait W to predict the individuals from the validation set. For that purpose, we extended the CDmean criterion [START_REF] Rincent | Maximizing the Reliability of Genomic Selection by Optimizing the Calibration Set of Reference Individuals: Comparison of Methods in Two Diverse Groups of Maize Inbreds (Zea mays L.)[END_REF] to a multi-trait CDmean criterion (CDmulti). We first computed the prediction error variance (PEV), using the following equations [START_REF] Lynch | Genetics and analysis of quantitative traits[END_REF]:

[ 𝑋 𝑇 (Σ ε -1 ⊗ 𝐼)𝑋 𝑋 𝑇 (Σ ε -1 ⊗ 𝐼)𝑍 𝑍 𝑇 (Σ ε -1 ⊗ 𝐼)𝑋 𝑍 𝑇 (Σ ε -1 ⊗ 𝐼)𝑍 + (Σ a -1 ⊗ 𝐾 -1 ) ] [ 𝛽 â] = [ 𝑋 𝑇 (Σ ε -1 ⊗ 𝐼)𝑦 𝑍 𝑇 (Σ ε -1 ⊗ 𝐼)𝑦 ] ( 9
)
where y is the vector of phenotypes, β is a vector of fixed effects (in our case it is the intercept), a is a vector of random genetic values, X and Z are the corresponding design matrices for the fixed and random effects respectively, K is the kinship matrix, Σa is the genetic variance-covariance matrix between traits, and Σε is the residual variance-covariance matrix between traits.

With the following notations:

[ 𝑋 𝑇 (Σ ε -1 ⊗ 𝐼)𝑋 𝑋 𝑇 (Σ ε -1 ⊗ 𝐼)𝑍 𝑍 𝑇 (Σ ε -1 ⊗ 𝐼)𝑋 𝑍 𝑇 (Σ ε -1 ⊗ 𝐼)𝑍 + (Σ a -1 ⊗ 𝐾 -1 ) ] -1 = [ 𝐶 11 𝐶 12 𝐶 21 𝐶 22 ] (10) we have 𝑃𝐸𝑉(â) = 𝑑𝑖𝑎𝑔(𝐶 22 ) (11)
The multi-trait CD was defined using the following equation:

𝐶𝐷 𝑚𝑢𝑙𝑡𝑖 (â) = 1 - 𝑃𝐸𝑉(â) 𝑑𝑖𝑎𝑔(Σ a ⊗𝐾) ( 12
)
𝐶𝐷 𝑚𝑢𝑙𝑡𝑖 is the expected reliability (before phenotyping) for each trait and each individual. It takes values between 0 and 1. The closer to 1 the criterion is, the more reliable the predictions are expected to be. We optimized the sampling by maximizing the mean of the 𝐶𝐷 𝑚𝑢𝑙𝑡𝑖 (𝐶𝐷 ̅̅̅̅ 𝑚𝑢𝑙𝑡𝑖 ) of individuals in the validation set for BMS, the trait of interest. For this a simple exchange algorithm was used. The initialization step consisted in randomly sampling calibration lines for BMS and/or W phenotyping and the corresponding 𝐶𝐷 ̅̅̅̅ 𝑚𝑢𝑙𝑡𝑖 was calculated for the validation set for BMS. For each next step, one of the two traits (BMS or W) was randomly chosen. Then a random exchange between one individual in the calibration set for the chosen trait and one individual which was not in the calibration set was realized. If this new calibration set resulted in a higher𝐶𝐷 ̅̅̅̅ 𝑚𝑢𝑙𝑡𝑖 , the new calibration set was accepted, otherwise the exchange was rejected. This procedure was repeated 2000 times to reach an optimum.

The genetic and residual variance-covariance matrices between traits (Σa and Σε) are estimated from phenotypic data. CDmulti is maximized before phenotyping. So those matrices have to be estimated using an independent dataset. To do so, we used 201 lines phenotyped for W and BMS by the French institute responsible for variety testing and registration (Groupement d'Etude et de contrôle des Variétés Et des Semences : GEVES) and genotyped using the same set of SNP markers.

We tested this optimization algorithm for several scenarios (Figure 1).

The optimization algorithm was performed with two different conditions:

1) in the first case, only individuals of the training set could be phenotyped for BMS and W (MT approach),

2) in the second case, only individuals from the training set could be phenotyped for BMS, while any of the individuals (from the training or the validation sets) could be phenotyped for the secondary trait W (TA approach).

In order to be able to compare GP models for different scenarios, lines from the validation set were never considered phenotyped for BMS, so that BMS predictive ability was always calculated on the same set of individuals.

We performed the optimization algorithm using different sample sizes that should be phenotyped for BMS and/or for W (Table 1).

We defined three budgets:

1) The most expensive budget allows to phenotype all the lines from the training set for BMS (47 700€),

2) The intermediate budget allows to phenotype all the lines from the training set for W and half of the lines for BMS (30 210€),

3) the less expensive budget allows to phenotype a quarter of the lines from the training set for BMS and three quarter of the lines from the training set for W (16 760€).

For a fixed budget, equation ( 8) has two unknowns (NBMS, NW), and there are many different possible combinations of individuals and traits to phenotype. We tested this optimization algorithm with 4 resource allocation strategies for each of the three budgets and under both MT and TA conditions. The resource allocation strategies correspond to a proportion (25, 50, 75 and 100%) of the training set (MT approach) or of the total population (TA approach)

that is phenotyped for W. Then the NBMS, the number of lines phenotyped for BMS, was calculated using the following equation obtained from equation ( 8):

𝑁𝑏 𝐵𝑀𝑆 = 𝐵-(𝐶𝑝 𝑊 ×𝑁𝑏 𝑊 ) 𝐶𝑝 𝐵𝑀𝑆 (13)
Training and validation sets were defined using 5-fold cross-validations. For each scenario (MT or TA approach, one budget and one resource allocation), we tested 20 folds (4 replicates x 5-fold cross-validations). In addition,

for each fold and each scenario the 𝐶𝐷 ̅̅̅̅ 𝑚𝑢𝑙𝑡𝑖 optimization algorithm was repeated 5 times, which means that for each scenario the procedure has been performed 100 times. The BMS BLUEs from the training set and W BLUEs from the training set (MT) and part of the validation set (TA) were used to predict BMS in the validation set. We compared predictive abilities obtained from optimized (𝐶𝐷 ̅̅̅̅ 𝑚𝑢𝑙𝑡𝑖 ) and random training sets (the mean predictability of 100 random sets for each scenario and each validation set).

Furthermore, the most expensive budget (47 700€) allowed to perform a ST approach in which all the individual from the training set were phenotyped for BMS. Therefore, we also compared predictive ability obtained for optimized training sets within each scenario with the predictive ability obtained with a ST model for BMS (calculated for each cross-validation that has been tested with the optimization algorithm).

In order to mimic real life situation, we also used a forward prediction strategy. In this case, the lines evaluated in 2003-2013 composed the training set and were used to predict the lines evaluated from 2014 to 2016.

The script of the optimization algorithm is available from the corresponding author on request. In addition, we extended the generalized CD of contrasts to the multi-trait context to adapt the optimization criterion to specific prediction objectives (see Appendix).

Results

Trait variation and variance components

Variance components and repeatabilities of yield, phenology and traits linked to bread making quality are presented in Table 2.

Plot repeatability ranged from 0.25 for crumb score to 0.81 for heading date. Design repeatability ranged from 0.66 for crumb score to 0.98 for heading date. Pearson's correlations between BLUEs of each trait are presented in Supplementary Figure S1. As expected, yield was negatively correlated with protein content (r = -0.63). The Chopin alveograph parameters were also correlated with each other. Indeed, Tenacity (P) was highly negatively correlated with Extensibility (L) (r = -0.6), while W and P were highly positively correlated (r = 0.63). Loaf volume, dough, crust and crumb scores are components of the BIPEA test which mechanically explains their correlation with BMS (the final BIPEA score).

As W is the trait that is the more correlated to BMS (r = 0.45) and its plot repeatability and design repeatability is higher (0.66 and 0.92 respectively) than for BMS (0.40 and 0.80 respectively), we used it as a secondary trait in MT and TA models in further analyses.

Contribution of glutenin (HMW-GS) markers to genomic prediction predictive ability

The 11 KASP markers derived from HMW-GS loci Glu-A1, Glu-B1, and Glu-D1 explained 31% and 8% of the phenotypic variation for W and BMS, respectively. The stepwise regression identified 3 markers explaining 30%

of the W genetic variance and 4 markers explaining 8% of the BMS genetic variance, thus nearly as much as the 11 markers used together.

In order to compare ST, ST-glu and ST-glu_stepwise models, only the 200 lines that were phenotyped and genotyped with the 35K array and the KASP markers derived from HMW-GS loci were used for genomic predictions (Figure 2). For ST model with 200 lines, predictive ability was 0.49 for W and 0.28 for BMS. ST-glu and ST-glu_stepwise models had almost the same predictive ability for both traits (0.55 for W and 0.28 for BMS).

They both improved the predictive ability of W but not that of BMS.

Single-trait versus multi-trait and trait-assisted genomic prediction models

Three models (ST, MT and TA) were tested using the 398 lines which were both genotyped and phenotyped for W and BMS. For ST model, predictive ability was 0.57 for W and 0.38 for BMS (Figure 2). Although MT model did not improve the quality of BMS predictions, TA model improved it (+0.15). In addition, we predicted W using BMS as correlated trait. MT model did not improve the quality of W predictions, unlike TA model. Besides, the gain in predictive ability generated by TA model was lower for W (+0.10) than for BMS (+0.15).

We compared MT and TA approaches with 3 different budgets and 4 different resource allocation strategies (Figure 3). TA method performed better in terms of predictive ability than MT method in each scenario in crossvalidation and in forward prediction. For each of the three budgets, the predictive ability for BMS was higher when all the lines (of the training set in MT case and of the total population in TA case) were phenotyped for the less expensive trait, W. The intermediate budget with at least 75 % of lines phenotyped for W was almost as accurate (0.5) than the expensive budget with 100% of lines phenotyped for W (0.51). This means that we can save 36% of the budget with this phenotyping strategy.

We used as a reference the predictive ability obtained with a ST model in which all the individuals from the training set were phenotyped for BMS (horizontal line in Figure 3). The cost of this ST approach is equal to the most expensive budget. For the intermediate budget, when all the lines from the training set were phenotyped for W, the predictive ability of BMS with MT model was close to the predictive ability obtained with ST model and the budget was reduced by 36%. In addition, for each budget, the predictive ability for BMS obtained using the TA model was higher than the predictive ability obtained with the ST model when at least three quarters of the lines were phenotyped for W. This means that TA approach improved the predictive ability compared to ST and MT approaches (maximum gain: +0.14 in cross-validation and +0.21 in forward prediction when all the lines were phenotyped for W) and allowed to reduce the cost by up to 65%.

In addition, the ST model to predict BMS for the lines evaluated between 2014 and 2016 using older lines (evaluated from 2003 to 2013) as training set, showed a lower predictive ability than in cross-validation scenarios (-0.2). MT approach slightly improved the predictive ability compared with ST model, but it remained under 0.22.

TA approach greatly improved this predictive ability (maximum predictive ability with random sampling: 0.41) whichever the budget. In addition, the predictive ability reached with TA approach in forward prediction was closed to the predictive ability obtained with ST approach in cross-validation scenario.

Selective phenotyping using Multi-trait CDmean criterion

To test whether we could minimize phenotyping or improve predictive ability by optimizing the choice of individuals to phenotype for BMS and/or W, we used a multi-trait CDmean criterion (𝐶𝐷 ̅̅̅̅ 𝑚𝑢𝑙𝑡𝑖 ) in an optimization algorithm. The multi-trait CD criterion allowed to improve slightly but systematically the predictive ability (+0.013 on average) for BMS compared to random sampling. The gain in terms of predictive ability allowed by the multitrait CD criterion was slightly higher for TA approach in forward prediction strategy (+0.023 on average).

On average, the gain in predictive ability generated by TA model compared to MT model was +0.13 in crossvalidation scenarios and +0.15 in forward prediction strategy. Note that for TA method using the multi-trait CDmean criterion, the number of lines phenotyped for W in the validation set was higher compared to random sampling (Supplementary Figure S2).

Impact of the cost ratio between W and BMS on BMS predictive ability

We evaluated the impact of the cost ratio between W and BMS on BMS predictive ability for MT and TA approaches. To do so, we compared predictive ability of BMS in several scenarios (Supplementary Table S2) by varying the cost ratio (ratio = 1/2, 1/3, 1/4, 1/5, 1/10 and 1/20). We fixed the budget allocated to phenotyping (the three budgets described above), the cost of phenotyping one line for BMS cost for one line (CpBMS=150€) and the percentage of the lines phenotyped for W (100% of the training set in MT approach and 100% of the total population in TA approach).

When the cost difference between the two traits was small and the budget allocated to phenotyping was low, the budget was not sufficient to phenotype each line for W (Supplementary Table S2). Therefore, only the results obtained with the intermediate and the expensive budget were shown in the Figure 4.

When the cost difference between the two traits increased, the number of lines phenotyped for BMS increased as well, and consequently the predictive ability of BMS was improved. Indeed for the expensive budget (which corresponds to the budget needed to phenotype each line of the training set for BMS), the predictive ability of BMS ranged from 0.35 when the ratio was 1/2 to 0.38 when the ratio was 1/20 for MT approach and ranged from 0.48 when the ratio was 1/2 to 0.52 when the ratio was 1/20 for TA approach. For the intermediate budget, the same trends were found. We could not evaluate the predictive ability of BMS with the TA model when the ratio between the costs of phenotyping of each trait was 1/2 because NBMS was too small and it led to convergence model issues. Even if the interest of using a TA approach rather than a ST approach was slightly lower when the cost difference decreased, TA model still performed better than ST approach when the ratio was 1/2 and for a fixed budget (gain in predictive ability: +0.09).

Impact of the number of phenotypic records on predictive ability

We evaluated the impact of the number of phenotypic records to estimate BMS and W BLUEs on BMS predictive ability for ST, MT and TA models (Figure 5). The number of phenotypic records for each line (except the control lines) ranged from 1 to 4. The three models were tested using the 322 lines which were both genotyped and phenotyped for W and for BMS in at least 4 environments.

We also evaluated the impact of control lines (evaluated 7 environments on average for BMS) to estimate BLUEs on BMS predictive ability. To do so, we compared two cases regarding controls: 1) all the phenotypic records of control lines were used to estimate BLUEs in order to have a more connected design; 2) we sampled the same number of control lines as candidate lines to calculate BLUEs.

For the three models, the predictive ability increased with the number of phenotypic records when we used all controls. The highest gain in terms of predictive ability was when the traits were evaluated in two environments compared to one (+ 0.11, + 0.08 + 0.07 for ST, MT, TA respectively). In addition, the predictive ability for using three records was slightly higher compared to two (+ 0.05 on average) and almost similar using four records (+ 0.02 on average) whichever the model used. For each scenario, the predictive ability was higher when all the phenotypic records of control lines were used to calculate BLUEs (+0.04, +0.04 +0.03 for ST, MT, TA respectively for 2 records). Note that the budget allocated to phenotyping can be significantly lower (-50 or 25% for 2 or 3 records). In addition, when both traits were measured in more than one environment, MT did not allow to improve predictive ability of BMS compared to ST model. But TA model always performed better than MT and ST models whichever the number of phenotypic records.

Discussion

This study focused on the comparison of different strategies to optimize resource allocation in breeding programs by using information from a cheap correlated trait to predict an expensive trait of interest. We used winter bread wheat baking quality (BMS) as a case study. Single-trait (ST), multi-trait (MT) and trait-associated (TA) strategies were compared in terms of predictive ability and cost. In multivariate approaches, information from dough strength (W) which is correlated (r = 0.45) to BMS and is less expensive to phenotype, was used. In MT analyses, only lines of the training set were phenotyped for W, whereas lines of the training and the validation sets were evaluated for W in TA analyses. In addition, we proposed and tested a multi-trait CDmean criterion (CDmulti) to optimize the choice of the set of lines to be phenotyped in the multi-trait context. To conduct this study, a public French breeding population of winter bread wheat lines was analyzed.

Ability of glutenin (HMW-GS) markers to predict BMS

In order to complement scarce phenotyping or pedigree information, [START_REF] Lande | Efficiency of marker-assisted selection in the improvement of quantitative traits[END_REF] proposed to rank individuals using molecular scores calculated as the sum of the effects of their alleles at selected QTLs. In this strategy, a threshold has to be chosen to determine the list of loci that significantly contribute to the trait variation. As it has been shown that allele effects were often overestimated in QTL detection [START_REF] Beavis | QTL analyses: power, precision, and accuracy[END_REF], and that the infinitesimal model (Fisher 1918) that considers many loci of small effects was the best model to explain the variation of many quantitative traits such as yield, [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF] proposed instead to use all available independent markers. To do so, phenotypic values are regressed in a linear model on all markers simultaneously, considering them as random effects with the same small genetic variance. The genetic values of individuals are estimated by best linear unbiased prediction. The critical difference with the Lande and Thompson approach is that this method does not require setting a significance threshold for the loci selected for trait prediction since all the markers are used in the model. In some cases however, explicitly modeling major QTLs as fixed effects in genomic prediction models can improve prediction accuracy in simulated data [START_REF] Bernardo | Genomewide Selection when Major Genes Are Known[END_REF]) and experimental bread wheat data [START_REF] Arruda | Comparing genomic selection and marker-assisted selection for Fusarium head blight resistance in wheat (Triticum aestivum L.)[END_REF][START_REF] Michel | Improving the baking quality of bread wheat by genomic selection in early generations[END_REF][START_REF] Sarinelli | Training population selection and use of fixed effects to optimize genomic predictions in a historical USA winter wheat panel[END_REF]. In bread wheat, HMW-GS are important determinants of rheological properties of dough (reviewed by [START_REF] Anjum | Wheat Gluten: High Molecular Weight Glutenin Subunits-Structure, Genetics, and Relation to Dough Elasticity[END_REF]). [START_REF] Ravel | SNP markers for early identification of high molecular weight glutenin subunits (HMW-GSs) in bread wheat[END_REF] developed a set of KASP markers derived from glutenin DNA sequences to identify HMW-GS alleles. They showed that these markers detected more alleles and explained more rheological traits variation than SDS-PAGE method, which was commonly used in plant breeding. In our study, we integrated 11 of these HMW-GS markers as fixed effects into ST models. It improved the predictive ability of W compared to ST model. These results are consistent with the conclusion of [START_REF] Michel | Improving the baking quality of bread wheat by genomic selection in early generations[END_REF] who found that prediction accuracy was increased for several dough rheological traits when Glu-1 loci markers were including as fixed effects into the RR-BLUP model. But this model did not improve the predictive ability of BMS. This result might be due to the low percentage of genetic variance of BMS explained by the HMW-GS markers (8%). [START_REF] Bernardo | Genomewide Selection when Major Genes Are Known[END_REF] actually demonstrated that only

QTLs that explained at least 10% of the genetic variance could improve prediction accuracy.

Prediction of unphenotyped individuals (MT)

Our results were consistent with some previous studies on experimental data in which MT models did not improve the predictive ability for individuals that were neither phenotyped for the trait of interest nor for the secondary trait, as in pine tree (Jia and Jannink 2012), soybean [START_REF] Bao | Association mapping and genomic prediction for resistance to sudden death syndrome in early maturing soybean germplasm[END_REF], rye [START_REF] Schulthess | Multiple-trait-and selection indices-genomic predictions for grain yield and protein content in rye for feeding purposes[END_REF]), maize (Dos Santos et al. 2016), wheat [START_REF] Schulthess | Advantages and limitations of multiple-trait genomic prediction for Fusarium head blight severity in hybrid wheat (Triticum aestivum L.)[END_REF][START_REF] Michel | Improving the baking quality of bread wheat by genomic selection in early generations[END_REF][START_REF] Lado | Resource allocation optimization with multi-trait genomic prediction for bread wheat (Triticum aestivum L.) baking quality[END_REF]) and sorghum [START_REF] Fernandes | Efficiency of multi-trait, indirect, and trait-assisted genomic selection for improvement of biomass sorghum[END_REF].

Previous studies showed that MT could be beneficial in terms of accuracy compared to ST model when a part of the individuals of the training set were not phenotyped for the targeted trait but the correlated trait was evaluated for a larger number of individuals [START_REF] Guo | Comparison of single-trait and multiple-trait genomic prediction models[END_REF][START_REF] Michel | Improving the baking quality of bread wheat by genomic selection in early generations[END_REF]). We found a similar trend in our study. In fact, when only half of the lines in the training set were phenotyped for BMS but all the lines (training and validation sets)were phenotyped for W, the predictive ability was similar compared to a ST model in which all the lines from the training set were phenotyped for BMS, but the budget was reduced by one third.

In the INRA-Agri-Obtentions bread wheat breeding program, F8 and F9 lines are evaluated in 9 environments in France on average. BMS is evaluated in 7 environments on average for control lines and in 4 environments for most of the other lines. Those measurements are very expensive (150€) and in this study, we tested if we could decrease the number of records (or the number of environments) without affecting the predictive ability of BMS.

As expected, predictive ability increased with the number of phenotypic records. The higher gap of predictive ability was observed when 2 phenotypic records instead of 1 were used to predict BMS with ST model (+0.08).

However, with 3 phenotypic records instead of 4, the predictive ability was little affected (+0.02), but the budget allocated to phenotyping was reduced by 25%. The predictive ability was higher (+0.03) when all the phenotypic records of control lines were used to calculate BLUEs. Controls are essential in the breeding program in order to take into account G×E and it is important to evaluate them in a large number of environments. By constrast, the number of records for other lines could be decreased with little effect on the predictive ability. Although phenotypic records were selected randomly for each line, we could further optimize the experimental design [START_REF] Heslot | Optimization of selective phenotyping and population design for genomic prediction[END_REF] and control lines distribution in particular.

Prediction when (part of) the candidates are phenotyped for the correlated trait (TA)

TA method consists of predicting a trait using a multi-trait model when at least some of the candidates are phenotyped for one or several correlated traits. Several studies using experimental data demonstrated that TA models perform better than ST and MT models in terms of accuracy. TA models using high-throughput phenotyping for instance improved the prediction accuracy of bread wheat grain yield by up to 70% [START_REF] Rutkoski | Canopy Temperature and Vegetation Indices from High-Throughput Phenotyping Improve Accuracy of Pedigree and Genomic Selection for Grain Yield in Wheat[END_REF][START_REF] Sun | Multitrait, Random Regression, or Simple Repeatability Model in High-Throughput Phenotyping Data Improve Genomic Prediction for Wheat Grain Yield[END_REF][START_REF] Crain | Combining High-Throughput Phenotyping and Genomic Information to Increase Prediction and Selection Accuracy in Wheat Breeding[END_REF]. TA models also improved bread wheat baking quality related parameters using protein content [START_REF] Michel | Improving the baking quality of bread wheat by genomic selection in early generations[END_REF] or dough rheological traits [START_REF] Lado | Resource allocation optimization with multi-trait genomic prediction for bread wheat (Triticum aestivum L.) baking quality[END_REF]) as correlated traits. Predictive ability of Fusarium head blight severity in hybrid bread wheat was improved using plant height and heading date as correlated traits [START_REF] Schulthess | Advantages and limitations of multiple-trait genomic prediction for Fusarium head blight severity in hybrid wheat (Triticum aestivum L.)[END_REF]. Wheat is not the only specie for which TA approach has been tested, [START_REF] Fernandes | Efficiency of multi-trait, indirect, and trait-assisted genomic selection for improvement of biomass sorghum[END_REF] showed that TA models increased prediction accuracy by up to 50% when using plant height as correlated trait to predict yield in sorghum.

Our results were consistent with these studies, since predictive ability using TA model was higher than ST and MT predictive abilities whichever the studied scenario.

For each studied budget (low, intermediate, high), the predictive ability of BMS was higher when all lines (training and validation) were phenotyped for the less expensive trait, W. It seems more important to phenotype a maximum of lines with a highly heritable trait than just a few lines for an expensive and less heritable trait. For similar budget, TA allowed higher predictive ability than ST (maximum gain: +0.14). For similar predictive ability (approximately 0.38), TA allowed reducing the budget by up to 65%. So TA is a way to improve genetic gain by cost unit, either by reducing phenotyping cost and budget or increasing the program's size for the same budget.

We showed that the gain in predictive ability increased with the difference of phenotyping cost between the expensive and the cheap trait. These results highlighted the advantage of using a correlated trait as cheap as possible in multi-variate genomic prediction models. Nevertheless, TA model was still interesting compared to ST model with the smallest cost difference (ratio = 1/2) tested in this study (gain in predictive ability: +0.09). In addition, even if it has no practical interest, we tried to predict W (plot repeatability = 0.66) using BMS (plot repeatability = 0.4) information and TA model, to compare the effect of the secondary trait's heritability (approximated by plot repeatability) on predictive ability. The gain was inferior to those when we predicted BMS using W. This is congruent with the hypothesis that the higher the heritability of the secondary trait compared to the trait of interest, the higher the gain in predictability [START_REF] Jia | Multiple-Trait Genomic Selection Methods Increase Genetic Value Prediction Accuracy[END_REF][START_REF] Hayashi | A Bayesian method and its variational approximation for prediction of genomic breeding values in multiple traits[END_REF][START_REF] Guo | Comparison of single-trait and multiple-trait genomic prediction models[END_REF].

Forward prediction

In order to study a situation more similar to the one the breeders are dealing with, we investigated a forward prediction strategy. ST predictive ability was almost divided by two when the oldest lines were used to predict the most recent lines compared to a cross-validation scenario. This discrepancy in predictive ability cannot be attributed to differences in the size of training sets. Indeed, the size of training set in cross-validation scenario was 318 lines vs 309 in the case of forward prediction. [START_REF] Storlie | Genomic Selection Accuracy using Historical Data Generated in a Wheat Breeding Program[END_REF] also reported that forward prediction in wheat is less accurate than cross-validation. Since the youngest lines are not differently (or less) genetically related to the oldest lines than the oldest lines with each other (results not shown), the most likely explanation is that new candidates are evaluated in new environments (years and locations) that have never been observed.

Although BLUEs were used to correct the environmental main effect, they did not correctly (efficiently) account for G×E interactions. Using TA by phenotyping some candidates with a correlated trait, we could take into account G×E in the new environment and better predict the candidates. This result suggests that TA model are particularly interesting to predict lines evaluated in new environments.

Moreover, a genomic Bayesian multi-trait and multi-environment model was developed by Montesinos-Lopez et al (2016) to take into account trait x genotype x environment interaction that could be interesting to test in that context.

Optimization of the training set based on a multi-trait CD criterion (CDmulti)

The size and the composition of the training population are key parameters of the accuracy of genomic prediction models. Optimizing the training set by selecting the most predictive individuals instead of using a random sampling leads to an increase of the predictive ability [START_REF] Rincent | Maximizing the Reliability of Genomic Selection by Optimizing the Calibration Set of Reference Individuals: Comparison of Methods in Two Diverse Groups of Maize Inbreds (Zea mays L.)[END_REF][START_REF] Akdemir | Optimization of genomic selection training populations with a genetic algorithm[END_REF]Isidro el al. 2015;[START_REF] Rincent | Predicting genomic selection efficiency to optimize calibration set and to assess prediction accuracy in highly structured populations[END_REF][START_REF] Sarinelli | Training population selection and use of fixed effects to optimize genomic predictions in a historical USA winter wheat panel[END_REF]. In this study, we specifically developed a multi-trait CDmean criterion (CDmulti), which slightly but systematically improved the predictive ability for BMS. The gain of predictive ability was a bit higher for TA approach than MT approach, in particular with the forward prediction. This may be due to the fact that 𝐶𝐷 ̅̅̅̅ 𝑚𝑢𝑙𝑡𝑖 leads to sample more lines to phenotype for W in the validation set because lines are more related to the candidates compared to the training set (Supplementary Figure S2). In this example the use of 𝐶𝐷 ̅̅̅̅ 𝑚𝑢𝑙𝑡𝑖 resulted in small increase of predictive ability, but we can suppose that it could be more valuable in datasets with different genetic structure, or in the situation of lower budget in which the choice of the phenotyped lines would be essential.

Note that the criterion CDmulti can be adapted to particular prediction objectives (eg: optimize predictive abilities within families as in [START_REF] Rincent | Predicting genomic selection efficiency to optimize calibration set and to assess prediction accuracy in highly structured populations[END_REF]) by using the corresponding contrasts (see Appendix for the generalized version of CDmulti).

Conclusion

In this study, we addressed the question of resource allocation, which is hardly ever mentioned in multi-trait studies. We demonstrated how to decrease the breeding cost of an expensive trait by using the information of a cheap correlated trait using trait-assisted (TA) genomic predictions. TA models are useful tools for breeders to optimize resource allocation while maintaining the predictive ability for the traits of interest. In our study, it reduced the budget allocated to phenotyping by up to 65% for Bread Making Score (BMS). TA models can also be used to improve the predictive ability of the trait of interest (+0.14 using cross-validation and +0.21 using forward predictions) for a fixed budget. The less expensive and labor-intensive, and the more heritable the correlated trait is compared to the trait of interest, the more TA could be useful in a breeding program.

In addition, we highlighted the importance of control lines that are phenotyped in a large number of environments and which allow to decrease the number of phenotypic records (3 instead of 4) for the other lines with a small effect on BLUEs estimation. MT_CDmulti, TA_CDmulti: MT and TA approaches using the CDmulti to sample lines to phenotype, respectively.

MT_Random, TA_Random: MT and TA approaches using a random sampling, respectively.

Appendix

We extended here the generalized CD [START_REF] Laloë | Precision and information in linear models of genetic evaluation[END_REF] to the multi-trait context. The objective is to compute the expected reliability (before phenotyping) in a multi-trait context for contrasts corresponding to the prediction objectives. To give an example, it is clear that if the objective is to accurately predict the difference (contrasts) between families, or to focus on predictions within families, the contrasts to be considered will be different and the optimal calibration sets as well. By defining contrasts corresponding to the prediction objectives, one can adapt the CDmulti criterion to the specific prediction objectives (see [START_REF] Rincent | Predicting genomic selection efficiency to optimize calibration set and to assess prediction accuracy in highly structured populations[END_REF] for more details in the single trait context).

For a given contrast c, the generalized Prediction Error Variance (PEV(c)) in the multi-trait context is equal to:

PEV(c) = 𝑐 𝑇 (𝑍 𝑇 𝑀𝑍 + (𝛴 𝑎 -1 ⊗ 𝐾 -1 )) -1 𝑐 (14)

With M = (𝛴 𝜀 -1 ⊗ 𝐼) -(𝛴 𝜀 -1 ⊗ 𝐼)𝑋(𝑋 𝑇 (𝛴 𝜀 -1 ⊗ 𝐼)𝑋) -1 𝑋 𝑇 (𝛴 𝜀 -1 ⊗ 𝐼)

where X and Z are the design matrices for the fixed and random effects respectively, K is the kinship matrix, Σa is the genetic variance-covariance matrix between traits, and Σε is the residual variance-covariance matrix between traits. ⊗ indicate the Kronecker product operator between matrices.

The generalized multi-trait CD for a given contrast c is equal to:

CDmulti (c) =

𝑐 𝑇 ((𝛴 𝑎 ⊗𝐾)-(𝑍 𝑇 𝑀𝑍+(𝛴 𝑎 -1 ⊗𝐾 -1 ))

-1

)𝑐 𝑐 𝑇 (𝛴 𝑎 ⊗𝐾)𝑐 [START_REF] Endelman | Ridge regression and other kernels for genomic selection with R package rrBLUP[END_REF] As a reminder, a contrast is a vector whose elements sum to 0 and indicating the difference in which we are interested. For example, if we are interested in accurately predicting the difference between individual 1 and individual 2, the contrast to consider will be: 𝑐 𝑇 = [1, -1,0,0,0,0, … ].

FigureFig. 1

 1 Figure

Fig. 2

 2 Fig.2 Comparison of predictive ability using single-trait and glutenin diagnostic markers (ST-glu and ST-glu_stepwise), single-trait (ST), multi-trait (MT) and trait assisted (TA) genomic prediction models N: Number of lines. On the left: comparison of predictive ability for each trait using three single-trait models (ST, ST-glu and ST-glu_stepwise). On the right: comparison of predictive ability using single-trait and multi-trait models. For the prediction of BMS, W was used as the correlated traits in bivariate models (MT and TA). For the prediction of W, BMS was used as the correlated traits in bivariate models (MT and TA). Error bars stand for standard deviation.

Fig. 3

 3 Fig. 3 Comparison of predictive ability for BMS in each scenario Results obtained with the 5-fold cross-validation method on the top. Results obtained when only the most recent lines are predicted on the bottom (Forward prediction). Solid lines correspond to predictive ability obtained using the ST model in cross-validation scenario. Dashed lines correspond to predictive ability obtained using the ST model in forward prediction scenario. Expensive, intermediate and cheap correspond to the three budgets.

Fig. 4

 4 Fig. 4 Impact of cost ratio between W (cheap) and BMS (expensive) traits on the predictive ability of BMS using multi-trait (MT) and trait assisted (TA) genomic prediction models The dashed line corresponds to predictive ability obtained using ST model. Ratio between phenotyping costs: 𝐶𝑝 𝑊 𝐶𝑝 𝐵𝑀𝑆 with 𝐶𝑝 𝑊 W cost for one line and 𝐶𝑝 𝐵𝑀𝑆 BMS cost for one line. CpBMS is a constant (150€) and CpW depends on the cost ratio. "Expensive" and "intermediate" correspond to the two budgets. The number on each barplot corresponds to the number of lines phenotyped for BMS. In each scenario, 100% of the lines of the training set (MT) or of the total population (TA) were phenotyped for W.

Fig. 5

 5 Fig. 5 Impact of the number of phenotypic records on predictive ability of BMS Case 1: all the phenotypic records of control lines were used to estimate BLUEs in order to have a more connected design. Case 2: phenotypic records of control lines were sampled in the same way as candidate lines to calculate BLUEs. In both cases, the total population (training and validation sets) includes 398 lines. Error bars stand for standard deviations.

phenotyped for W (In the training set for MT and in the total population for TA)

  The tools and approaches developed in this study and exemplified using BMS prediction of bread wheat can be used in any species when relevant cheap correlated traits are available.

			MT		TA
	Budget	% of lines NW	NBMS	NW	NBMS
		100	318	275	398	264
	47 700€	75	238	286	298	278
		50	159	296	199	291
		25	79	307	99	304
		100	318	159	398	148
	30 210€	75 50	238 159	169 180	298 199	161 174
		25	79	190	99	188
		100	318	69	398	58
	16 760€	75 50	238 159	80 90	298 199	72 85
		25	79	100	99	98

The multi-trait CD criterion (CDmulti) specifically developed in this study allowed a slight but systematic improvement of BMS predictive ability. It would be interesting to evaluate this criterion in a larger and more diverse population.

Table 1

 1 Number of selected individuals in each scenario NW, NBMS: Number of lines phenotyped for W and BMS, respectively

			Trait	Mean	Var	𝝈 𝑬 𝟐	𝝈 𝑮 𝟐	𝝈 𝑮𝒙𝑬 𝟐	𝟐 𝝈 𝒆	Plot repeatability	Design repeatability
		Grain yield	91.34	269.67	236.06	14.05	26.55	9.27	0.28	0.84
		Plant height	87.85	94.27	50.85	30.31	5.87	4.94	0.74	0.95
		Heading date	140.98	63.04	45.47	8.08	7,72x10 -4	1.93	0.81	0.98
		Protein content	11.17	1.89	1.4	0.26	0.21	0.1	0.46	0.91
	Chopin alveograph	Parameters	Dough strength:W Tenacity: P Extensibility: L 85.05 205.63 3844.01 77.46 695.31 936.17 𝑷 𝑳 ⁄ ratio 1.16 0.74	788.98 2030.03 166.1 413.3 254.77 465.46 0.14 0.34	------------	1045.85 117.67 282.5 0.24	0.66 0.78 0.62 0.59	0.92 0.95 0.9 0.89
			Dough score	79.24	215.21	27.78	69.5	---	120.13	0.37	0.77
	Test	Crust score	50.18	483.88	96.94	145.1	---	249.32	0.37	0.77
	BIPEA	Crumb score Loaf volume 1482.44 48629.99 10647.39 15815.57 93.98 61.16 15.13 12.16	------	36.53 23975.46	0.25 0.4	0.66 0.79
			BMS	223.4	1276.6	141.13	466.23	---	688.96	0.4	0.8

Table 2

 2 Summary statistics, variance components and repeatabilities for the main traits and traits linked to bread

	making quality
	σ 2 G; σ 2 G×E and σ 2 ε are the genotype, genotype x environment and residual variances respectively

Acknowledgements

The authors thank the genotyping platform GENTYANE at INRA Clermont-Ferrand (gentyane.clermont.inra.fr) which has conducted the genotyping. The work in experimental units by INRA (Clermont-Ferrand, Dijon, Estrées-Mons, Lusignan, Le Moulon and Rennes) and in Agri-Obtentions is also gratefully acknowledged. Doctoral work of SBS was funded by a grant from the Auvergne-Rhônes-Alpes region and from the European Regional Development Fund (FEDER).

This work was supported by the Breedwheat project thanks to the funding from the French Government managed by the National Research Agency (ANR) in the framework of Investments for the Future (ANR-10-BTBR-03)

France AgriMer and the French Fund to support Plant Breeding (FSOV).

Compliance with ethical standards

Author contribution statement

JA, FXO, BR and EH designed the field trials and collected the phenotypic data. RR supported in statistical analysis and in developing the multi-trait CDmean algorithm. CR developed KASP markers derived from HMW-GS loci. SBS analyzed the data and wrote the manuscript. SB and GG guided through the study and helped improving the manuscript. All authors approved the final manuscript.

Conflict of interest

The authors declare no conflict of interest

List of tables

Supplementary material

Name of trial GPS coordinates

Champagne-Céréales


E Kume

L Noirez

From Linear to Harmonic Thermal response of mesoscopic liquids upon mechanical shear strain

In the conventional picture and at the macroscopic scale, the temperature of a liquid bath in the quiescent state is uniform down to density fluctuation length scales. Here we examine the impact of a low frequency shear mechanical field (Hz) on the thermal equilibrium of mesoscopic liquids (Polypropylene glycol and pentadecane away from any phase transition). We show the emergence of both cooling and heating shear waves of several tens microns widths synchronous with the applied shear strain wave. The thermal wave is stable at relatively low strain amplitude and low frequency while thermal harmonics develop by slightly increasing the frequency or the strain amplitude. The liquid layer behaves as a thermoelastic medium. These results seem to confirm recent theoretical models predicting that liquids support shear elastic waves up to finite propagation length scale of the order the thermal wave.

Introduction:

The mechanical properties of fluids are extensively studied in terms of stress response but not in terms of thermal response. Mechanically induced thermal effects are typically solid-like properties, known as thermoelasticity. For small changes of temperatures and volumes, the thermoelastic effects are generally weak, reversible and giving rise to an adiabatic transformation where the compressibility/dilatation ability of the material is varied under mechanical stress: T  -

 𝐿 𝐶 
.  where  is the stress component, L the linear thermal expansion and C the specific heat at constant stress. This thermal variation can be positive or negative depending on the nature of the transformation leading to stretched or compressed regions respectively [START_REF] Biot | Thermoelasticity and Irreversible Thermodynamics[END_REF]. Without heat source, no temperature variation is predicted for fluids upon shear mechanical action as the energy is expected to be entirely dissipated on large volumes [START_REF] Thomson | Kinetic Theory of the Dissipation of Energy[END_REF]. The situation might be different by exploring low scale liquid behaviour. Indeed while the liquid definition implies an absence of static shear elasticity, shear elasticity has been recently predicted and experimentally measured in mesoscale liquids [START_REF] Derjaguin | The complex shear modulus of polymeric and small-molecule liquids[END_REF][START_REF] Noirez | Identification of finite shear-elasticity in the liquid state of molecular and polymeric glass-formers[END_REF][START_REF] Noirez | Revealing the solid-like nature of glycerol at ambient temperature[END_REF][START_REF] Noirez | Identification of a low-frequency elastic behaviour in liquid water[END_REF][START_REF] Chushkin | Low-frequency elastic behavior of a supercooled liquid[END_REF]. Volino proposed a non-extensive model to describe dimensiondependent viscoelasticity of fluids [START_REF] Volino | Théorie visco-élastique non-extensive[END_REF]. Revisiting the Frenkel model, Trachenko et al. introduce the concept of local dynamic compressive stress and revisits the Maxwell approach showing that a solidlike approach can be treated on equal footing introducing the notion of finite shear wave propagation length (gapped momentum states) [START_REF] Trachenko | Lagrangian formulation and symmetrical description of liquid dynamics[END_REF][START_REF] Trachenko | Quantum dissipation in a scalar field theory with gapped momentum states[END_REF][START_REF] Baggioli | Gapped momentum states[END_REF][START_REF] Zaccone | Explaining the low-frequency shear elasticity of confined liquids[END_REF][START_REF] Zaccone | Approximate analytical description of the nonafne response of amorphous solids[END_REF]. These models converge in interpreting the macroscopic liquid behavior as an asymptotic branch where, at a smaller scale, liquid molecules behave elastically with suppressed relaxations and where the dynamic role of intermolecular interactions is central [START_REF] Zaccone | Explaining the low-frequency shear elasticity of confined liquids[END_REF][START_REF] Zaccone | Approximate analytical description of the nonafne response of amorphous solids[END_REF]. Since evidence of small-scale elasticity seems established, it is natural to test whether another solid property, thermoelasticity, is also observable. Thermal variation involves an energy change between molecules through short time dependent interactions of the vibrations of the molecules. The far-infrared portion of the wide electromagnetic spectrum covers wavelengths from approximately 3 to 14 µm. This range concerns typically the intermolecular vibrations of the liquid dynamics; i.e. it characterizes the way molecules vibrate in the "quasi-cage structure" of the neighbored molecules. Except for solid stress analysis [START_REF] Greene | Thermoelastic Stress Analysis[END_REF], this fundamental property that is the thermal radiation naturally emitted by the condensed matter above zero degree Kelvin is not exploited to access the liquid thermodynamic state. We present a high accuracy real-time micro-thermal analysis of a liquid layer of a few hundred of microns set under oscillatory shear deformation conditions. The studied liquid is the glass former polypropylene glycol (PPG) at room temperature; i.e. far from any phase transition. This liquid is nearly transparent to infrared wavelength so that the bulk behavior can be probed. We show that the liquid responds thermally to the low frequency mechanical excitation challenging the assumption of fast thermal relaxation in liquids. We characterize the thermal signal and identify thermal harmonics when the excitation becomes larger in magnitude and frequency. Preliminary tests were also carried out on single molecule non polar solvent, the pentadecane (alkane), allowing us to pretend that the following results might be representative of a generic liquid property.

Methods:

For the experiments, two ordinary liquids, the polypropylene glycol (PPG, Sigma-Aldrich Mn = 4000) which is a H-bond glass former with an extremely low evaporation rate, and the pentadecane (Sigma-Aldrich) which a molecular van der Waals liquid were used. The measurements took place at room temperature ~20°C, which is far from the phase transition of PPG (Tg = 180 K) and pentadecane (melting temperature: 282.3 K [START_REF] Trachenko | Quantum dissipation in a scalar field theory with gapped momentum states[END_REF]).

The liquids were confined in between two parallel high-energy a-Alumina plates. The plates were cleaned by heating at 400°C for an optimal liquid/substrate wetting ensuring ideal conditions for the transmission of a dynamic shear strain generated from the sinusoidal movement of the bottom plate  = 0 . sin(.t) where 0 is the strain amplitude and  the frequency. Imposing a sine wave deformation allows by changing the frequencies and amplitudes to examine the stability of the signal response. A strain-controlled rheometer (TA-Instruments ARES) was used to impose the frequency, the strain amplitude and the gap thickness a Keithley multimeter is coupled as analog output to the rheometer and to a PC to allow the synchronization with an infrared microbolometer camera. The camera matrix of 388 * 282 pixels of 20 µm uses wavelengths from 8 to 14 µm at a frequency reading of 27 Hz, making possible the observation of the thermal response during a full mechanical oscillation since the applied frequency varies from 0.5 to 5 rad/s (0.08 to 0.8 Hz). The measurements were conducted with a thermal sensor placed at 50mm away from the plate-liquid-plate plane focusing at a depth of ~0.85mm inside the liquid, making possible to access the bulk of the liquid. The infrared data were corrected by subtracting the average thermal frame recorded at the equilibrium from every frame. The kinetic representation of the thermal mapping was created by the compacted succession of each frame next to each other.

Results

General features:

Fig. 1a andb show a typical thermal evolution along two oscillatory periods of the liquid PPG at  = 0.5 rad/s and 1 rad/s respectively (at 240µm gap thickness and 0 =4000 %). This kinetic representation of the temperature in the gap reveals a succession of instant thermal states at different stages of the oscillatory deformation. The different colors in Fig. 1a and Fig. 1b indicate that the liquid undergoes a thermal change exhibiting both local negative and positive temperature variations that reach sizeable amplitudes of about  0.2° C. These thermal waves are reproducible, reversible and stable in time as long as the mechanical strain wave is applied. The coexisting hot and cold waves are roughly distributed in three bands (of ~ 50μm thickness each) that we name "bottom, middle and upper bands". The bands exhibit in phase and opposite thermal waves which compensate for each other; Fig. 1c and1d illustrates the temperature variation in each chosen bands. Cold and warm zones alternate synchronously with the applied excitation; i.e. while the central band is cooling down, the neighbored ones are heating up and reciprocally. At low strain values (< 500%), the signal of thermal response is weak and noisy since the variation is close to experimental error bar (~ 0.02°C), making difficult to characterize in detail the thermal response at low strain amplitude.

Figure 1. By applying a low frequency mechanical stimulus (~ Hz), the liquid emits a modulated thermal signal synchronous with the stimulus (here PPG-4000 confined in a 240µm gap (gap view), alumina substrate, room temperature measurements). a) Real-time mapping of the temperature γ = 4000% (gathering about ~800 frames): at ω = 0.5 rad/s b) same at ω = 1 rad/s. The dotted line is an eye guide for the applied strain. c) Thermal waves recorded for the bottom (), middle () and upper () bands respectively at γ = 4000% and  = 0.5 rad/s. The dotted line illustrates the applied shear strain. φ is the phase shift between the thermal and strain wave (φ  45° at  = 4000%). d) Same for  = 1 rad/s. e) Strain dependence of the maximum of the temperature variation T0(°C) versus shear strain and inserts representing the "gain" T/ versus shear strain for the bottom "hot" band: (•), the middle "cold" band: (■) and the upper band: (♦) respectively at ω = 0.5 rad/s, (f) same at  = 1 rad/s.

Strain dependence of the thermal signal:

The strain dependence of the cold wave is shown in Fig. 1e and Fig. 1f at 0.5 rad/s and 1 rad/s respectively. At 0.5 rad/s, the strain dependence of the main waves (hot and cold) is relatively simple exhibiting a nearly linear dependence over a large strain range (Fig. 1e). The thermal modulation can be modelled by a sin wave: T = T0 . sin(.t + δφ) where δφ is the phase shift that is nearly in phase at low strain ( < 1000% not shown), and above, reaches approximately /4. The sin wave is kept up to the highest reachable strain values (Fig. 1c is obtained at 4000%), indicating the additional energy associated to the increase of the shear strain is transferred without distortion to amplify the thermal wave. The amplitude of the thermal wave can be modeled by a linear strain-dependence reaching a maximum variation of ~  0.2°C (10 times over the accuracy) (Fig 1e). The linear relationship suggests that the thermal mechanism is occurring as soon as the smallest shear strain values (measurements at  < 100% are below the accuracy); i.e. inherent to the liquid deformation. The insert displays the rate of the temperature variation over the strain. ΔΤ/γ is equivalent to the "gain" which measures the quality of the energy conversion. This "gain" is nearly constant for the three bands (insert of Fig. 1e). Another remarkable feature is the symmetry of the evolution of the cold and hot bands pointing out a rigorous compensation to keep the net temperature unchanged.

At 1 rad/s, the strain dependence of the thermal signal is more complex. Two thermal regimes are observed depending on the strain amplitude:

-At relatively "low" strain amplitudes ( < 2500% at  = 1 rad/s), the thermal modulations obeys a sin wave with a phase shift smaller than /4 in advance with the applied strain wave (Fig. 2a). Specifically, the phase shift at  = 1000% is δφmiddle ~ 32° ± 2°, δφbottom ~ 29° ± 2° for the cold and the hot bands respectively. The liquid exhibits an output thermal response synchronous with the input strain wave. The "hot" and "cold" waves evolves symmetrically confirming that the system works nearly in an exact thermal compensation as the strain amplitude increases.

-At higher strain amplitude of the strain wave ( > 2500% at  = 1 rad/s), the linear relationship with the input strain wave does not hold anymore, meaning that non-linear effects start to take place. Fig. 1f shows the progression of the amplitude of both hot and cold wave amplitudes. For the highest strain values, the thermal variation remains constant as shear strain increases (high strain plateau). Simultaneously the shape of the thermal wave evolves and deforms. Fig. 2 details the evolution of the cold wave. The wave splits in harmonics, while increasing the thermal amplitude for higher strain values. The extrapolation to infinite strain amplitudes is a contraction of the thermal amplitude and possibly a vanishing of the thermal variation (experimental restrictions do not allow us to go further). The infinite strain corresponds to steady-state conditions; i.e. to the flow regime. The vanishing of the thermal behavior at infinite strain is in agreement with the common observation of absence of noticeable thermal effects at macroscopic scale.

Frequency dependence of the thermal signal:

Fig. 3a-3c shows the influence of the frequency on the thermal wave in a range of 0.5 -5 rad/s (measured at the maximum strain wave amplitude  = 4000%). At low frequency, the thermal waves reproduce the input strain waveform (Fig. 3a illustrates the cold band at ω = 0.5rad/s). As the frequency is increased, the thermal wave progressively loses the input sinus shape (Fig. 3b and3c) and its amplitude decreases of nearly 40% from 0.5 to 5 rad/s. The thermal signal was decomposed with a process known as Fast Fourier Transformation (FFT) [START_REF] Baggioli | Gapped momentum states[END_REF] to compare its frequency domain with the one of the applied shear strain signal. FFT analysis (Fig. 3d-3f) shows that higher order harmonics appear at high frequencies in addition to the fundamental frequency (Fig. 3b 3d-f). The thermal waves were fitted using the sum of the harmonics: 𝑓(𝑡) = 𝐴 0 + ∑ 𝛢 𝑖 sin (𝜔 𝑖 𝑡 + 𝜑 𝑖 ) 𝑛 𝑖=1

,where i is a numerator and n represents the number of harmonics used, Ai constant and φi the phase shift. The generation of thermal harmonics implies that the excitation frequency is faster than the thermal relaxation between two successive oscillations and so, the thermal response is no more stable.

3.

Influence of the frequency on the evolution of the cold thermal wave at large strain amplitude (γ = 4000%) and corresponding FFT analysis. Increasing the frequency produces thermal harmonics and a collapse of the thermal variation. Data points (blue circles) corresponds to the thermal variation ΔΤ within PPG-4000 middle band at 0. 240mm,,a

Probing van der Waals liquid:

Figure 4 demonstrates that the thermal response is also measurable for a very low molecular weight liquid; the pentadecane. We observe a weak thermal response for low gap thickness and high shear strain values with the same characteristics as the PPG, namely the generation of hot and cold bands. Specifically, thermal variation is not observed in the whole liquid, but near the walls (30 -40 μm width). The middle part of pentadecane gap, which is less effected than the surfaces, does not exhibit thermal oscillations. The weaker thermal response could be attributed to the non-polar nature of pentadecane, which is a van der Waals molecule while PPG is a hydrogen bonded oligomer. Previous studies report on liquid shear elasticity is higher and more easily identified in polymer melts [START_REF] Noirez | Identification of finite shear-elasticity in the liquid state of molecular and polymeric glass-formers[END_REF] than less viscous liquids (glycerol, water) [START_REF] Noirez | Revealing the solid-like nature of glycerol at ambient temperature[END_REF][START_REF] Noirez | Identification of a low-frequency elastic behaviour in liquid water[END_REF], implying that different cohesion of the liquids matter on the intensity of solid-state properties in liquids. Another possible cause of different magnitude of the thermal effect could be that the liquids have different natural frequencies for the thermal effect to manifest greater. On the other hand, the identification of thermal response in a non-polymeric liquid suggests that the observed thermo-elastic effect is generic to liquids.

Discussions-conclusions:

The present results validate the importance of a thermal approach for the identification of hidden liquid properties. We observe the emission of a synchronous hot and cold thermal modulation upon applying an oscillatory shear strain wave. The examination of the thermal wave versus frequencies and strain amplitude showed that the thermal response is a stable sin wave within a low frequency range ( 0.5 rad/s) at small gap (250µm). A nearly linear variation of the temperature variation is observable for both cold and hot thermal bands at moderate shear strain amplitude (Fig. 1e-1f). At high frequencies or at very large shear strain amplitudes, the thermal response exhibits a distortion of the wave indicating that the thermal signal is no more stable. The thermal waves exhibit a phase shift of about /4, a part of the shear energy is dissipated in the liquid. For frequency larger than 0.5 rad/s, the generation of thermal harmonics suggests that further dissipation takes place due secondary effects and an applied excitation faster than the relaxation time of the thermal effect (thus estimated at   1-2s). We identify odd and even harmonics for the thermal response as the strain and/or frequency increases. The extrapolation at high frequency or at high strain amplitude indicates a smoothing of the thermal oscillation towards a flat signal (constant temperature).

The thermal behavior reveals that in apparent "viscous" liquids, a more complex mechanism governs the mesoscale. It reveals the liquid ability to convert the shear wave energy in local thermodynamic states. Such effect has characteristics of a conservative mechanism: -the linearity as the shear strain increases, -synchronism and nearly instant thermal responses to the strain wave, and at the probed frequencies, an absence of heat conduction between bands and with the environment. These properties might be regarded as mechanically induced adiabatic thermodynamic events (at least in the linear domain). These are typical of a thermoelastic effect; i.e. mediated by a liquid elasticity varying as a function of temperature. Shear elasticity was indeed identified in mesoscale liquids by stress measurements and depending on the scale considered. Since the early works of Derjaguin in 1989 [START_REF] Derjaguin | The complex shear modulus of polymeric and small-molecule liquids[END_REF], evidence of the existence of long range interactions in liquids accumulates from both experimental results [START_REF] Derjaguin | The complex shear modulus of polymeric and small-molecule liquids[END_REF][START_REF] Noirez | Identification of finite shear-elasticity in the liquid state of molecular and polymeric glass-formers[END_REF][START_REF] Noirez | Revealing the solid-like nature of glycerol at ambient temperature[END_REF][START_REF] Noirez | Identification of a low-frequency elastic behaviour in liquid water[END_REF][START_REF] Chushkin | Low-frequency elastic behavior of a supercooled liquid[END_REF][START_REF] Shelton | Long-range orientation correlation in dipolar liquids probed by hyper-Rayleigh scattering[END_REF][START_REF] Lv | Viscosity of water and hydrocarbon changes with micro-crevice thickness[END_REF][START_REF] Li | Nonlinear Viscoelastic Dynamics of Nanoconfined Wetting Liquids[END_REF] and fundamental models. The later ones come from new theoretical developments [START_REF] Trachenko | Lagrangian formulation and symmetrical description of liquid dynamics[END_REF][START_REF] Trachenko | Quantum dissipation in a scalar field theory with gapped momentum states[END_REF][START_REF] Baggioli | Gapped momentum states[END_REF][START_REF] Zaccone | Explaining the low-frequency shear elasticity of confined liquids[END_REF][START_REF] Zaccone | Approximate analytical description of the nonafne response of amorphous solids[END_REF], which revisit the Frenkel theory of propagation of transverse waves in liquids [START_REF] Frenkel | Kinetic theory of liquids[END_REF]. It is also predicted that the low frequency shear modulus scales as G′ ~ h -3 [START_REF] Zaccone | Explaining the low-frequency shear elasticity of confined liquids[END_REF], explaining why elastic effects become dominant as the scale decreases. The present thermal results are an indirect proof of shear waves propagation in a liquid medium and shows that the accumulated shear energy modifies the liquid thermodynamic at small scale (<500μm). From a mechanical point of view, reversible temperature changes induced by the deformation known as thermo-(shear)elastic effect are forbidden by the conventional liquid incompressibility assumption. Non-zero thermoelastic coefficients (here estimated at Kshear= ΔT/γ about 3.6 ± 0.3 .10 -3 K -1 ) imply that isobaric expansion and isothermal compressibility of liquids are also non-zero, opening the route to very interesting and novel types of developments in the emerging fields of nano and microfluidics [paper in progress].

Figure 2 .

 2 Figure 2. Evolution of the cold thermal wave versus shear strain at constant frequency (1 rad/s) for 1000%, 2500% and 4000% shear strain amplitude (cyan, blue and dark blue points respectively). At large strain, the thermal wave exhibits harmonics with a smoothed amplitude. Orange dotted lines are fits following eq. 𝑓(𝑡) = 𝐴 0 + ∑ 𝛢 𝑖 sin (𝜔 𝑖 𝑡 + 𝜑 𝑖 ) 2 𝑖=1

  ) ω = 0.5 rad/s, b) ω = 1 rad/s, c) ω = 5 rad/s. Orange lines are fitting of eq. 𝑓(𝑡) = 𝐴 0 + ∑ 𝛢 𝑖 sin (𝜔 𝑖 𝑡 + 𝜑 𝑖 ) 𝑛 𝑖=1 for a) n=1 b) n=2 c) n=3, Fourier Transform signal of the waves depicted of d) Fig. 3a, e) Fig. 3b, f) Fig. 3c. Inset of each graph corresponds to the equivalent 2D thermal mapping.

Figure 4 .

 4 Figure 4. Thermal response of pentadecane (0.145mm, 1rad/s, 5000%). Red and blue points represent top and bottom band respectively. Black dotted line represents the applied strainalumina substrate, room temperature measurements.

Acknowledgments: We thank P. Baroni for instrumental innovation and assistance. This work has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 766007 and from the LabeX Palm (ANR-11-Idex-0003-02).