Paulino Pérez-Rodríguez

José Crossa
email: j.crossa@cgiar.org

Jessica Rutkoski

Jesse Poland

Ravi Singh

Andrés Legarra

Enrique Autrique

Gustavo De Los Campos

Juan Burgueño

Susanne Dreisigacker

SINGLE-STEP GENOMIC AND PEDIGREE GENOTYPE × ENVIRONMENT INTERACTION MODELS FOR PREDICTING WHEAT LINES IN INTERNATIONAL ENVIRONMENTS

Keywords: wheat pedigree prediction, wheat genomic prediction, single-step prediction accuracy, genotype × environment interaction, international environments

Genomic prediction models have been commonly employed in plant breeding but only in reduced data sets comprising a few hundred genotyped individual plants. However, pedigree information for an entire breeding population is frequently available, as are historical data on the performance of a large number of selection candidates. The single-step method extends the

genomic relationship information of genotyped individuals to pedigree information from a larger number of phenotyped individuals so as to combine relationship information on all members of the breeding population. Furthermore, genomic prediction models that incorporate genotype × environment interaction (G×E) have produced substantial increases in prediction accuracy compared to single-environment genomic prediction models. The main objective of this study was to show how to use single-step genomic and pedigree models to assess the prediction accuracy of a large number of CIMMYT wheat lines (58,798) evaluated in several simulated

INTRODUCTION

Global wheat production is increasing less than 1% annually, and recently wheat yields have stagnated in many regions of South Asia [START_REF] Ray | Recent patterns of crop yield growth and stagnation[END_REF]. In South Asia, the wheat crop is already being grown under high temperature conditions; however, due to climate change, temperatures could increase well beyond the optimal for growing wheat, which would further reduce grain yield. As a result, South Asian countries may not be able to meet the region's already growing demand for wheat grain.

Well managed crop improvement programs are necessary to increase food production in different parts of the world. Several molecular marker methods have proven their relevance in different cereal crops. Genomic selection (GS) is becoming a standard approach to achieve genetic progress in plants because it reduces the generation interval by reducing the need to have progeny field-tested every cycle. Breeding values can be predicted as the sum of the effects of all markers by regressing the values of the phenotypes on all markers [START_REF] Meuwissen | Prediction of total genetic values using genome-wide dense marker maps[END_REF].

Several authors have successfully implemented GS in plant breeding with intermediate-to-high density marker coverage for traits such as grain yield, biomass yield, resistance to several diseases, and flowering evaluated under different environmental conditions. Studies have demonstrated that some of the factors determining prediction accuracy in GS are the heritability of the trait, the number of markers, the size of the training population, the relationship between Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

3 the training and the testing sets, and genotype × environment interaction (G×E) [START_REF] De Los Campos | Predicting quantitative traits with regression models for dense molecular markers and pedigree[END_REF][START_REF] Crossa | Prediction of genetic values of quantitative traits in plant breeding using pedigree and molecular markers[END_REF][START_REF] Crossa | Genomic selection and prediction in plant breeding[END_REF][START_REF] Pérez-Rodríguez | Comparison between linear and non-parametric models for genomeenabled prediction in wheat. G3: Genes, Genomes[END_REF][START_REF] Burgueño | Genomic prediction of breeding values when modeling Genotype × Environment interaction using pedigree and dense molecular markers[END_REF][START_REF] Hickey | Factors affecting the accuracy of genotype imputation in populations from several maize breeding programs[END_REF][START_REF] González-Camacho | Genome-enabled prediction of genetic values using radial basis function[END_REF][START_REF] Riedelsheimer | Genomic and metabolic prediction of complex heterotic traits in hybrid maize[END_REF][START_REF] Weber | Efficiency of managed-stress screening of elite maize hybrids under drought and low nitrogen for yield under rainfed conditions in Southern Africa[END_REF]. Furthermore, including high-density marker platforms with G×E interactions adds power to GS models [START_REF] Burgueño | Genomic prediction of breeding values when modeling Genotype × Environment interaction using pedigree and dense molecular markers[END_REF]Jarquín et al., 2014;[START_REF] López-Cruz | Increased prediction accuracy in wheat breeding trials using a marker × environment interaction genomic selection model[END_REF][START_REF] Heslot | Genomic selection in plant breeding: a comparison of models[END_REF].

Recently, genomic predictions have been extensively studied in bread wheat using elite germplasm sets (de los [START_REF] De Los Campos | Predicting quantitative traits with regression models for dense molecular markers and pedigree[END_REF]Campos et al., , 2010;;[START_REF] Crossa | Prediction of genetic values of quantitative traits in plant breeding using pedigree and molecular markers[END_REF]Gónzalez-Camacho et al., 2012;[START_REF] Heslot | Genomic selection in plant breeding: a comparison of models[END_REF][START_REF] Pérez-Rodríguez | Comparison between linear and non-parametric models for genomeenabled prediction in wheat. G3: Genes, Genomes[END_REF][START_REF] López-Cruz | Increased prediction accuracy in wheat breeding trials using a marker × environment interaction genomic selection model[END_REF]. The results have proven that the use of dense molecular markers coupled with pedigree information increases the prediction accuracy of unobserved phenotypes. One of the problems usually encountered by GS in animal and plant breeding is that the number of evaluated lines exceeds the number of genotyped lines, due to genotypic costs. [START_REF] Nejati-Javaremi | Effect of total allelic relationship on accuracy of evaluation and response to selection[END_REF] were the first to propose incorporating genotypic information for predicting the breeding values of animals in a manner similar to the way pedigree information is used in the Best Linear Unbiased Predictor (BLUP) method. When the pedigrees of all phenotyped individuals were available but only some were genotyped, dairy cattle researchers [START_REF] Misztal | Computing procedures for genetic evaluation including phenotypic, full pedigree, and genomic information[END_REF][START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF][START_REF] Aguilar | Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF][START_REF] Aguilar | Efficient computations of the genomic relationship matrix and other matrices used in single-step evaluation[END_REF][START_REF] Christensen | Single-step methods for genomic evaluation in pigs[END_REF] derived a unified (single-step) computation approach for Genomic Best Linear Unbiased Predictor (ssGBLUP) for combining phenotypic, pedigree and genomic information based on [START_REF] Henderson | Rapid method for computing the Inverse of a Relationship Matrix[END_REF][START_REF] Henderson | A simple method for computing the inverse of the numerator relationship matrix used in prediction of breeding values[END_REF] standard mixed model equations.

These authors augmented pedigree-based relationship matrix (A) with contributions from genomic relationship matrix (G) of the genotyped individuals. They showed how to modify the original A matrix to obtain an H matrix that includes not only the pedigree-based relationship matrix but also a matrix that contains the differences between genomic-based and pedigree-based matrices. These authors also developed efficient computer algorithms for inverting matrix H computed from large numbers (millions) of animals in the data.

Although augmenting matrix A by using only a fraction of the individuals that were genotyped would reduce genotyping costs, the ssGBLUP method has not been extensively applied in plant breeding. Just recently, [START_REF] Ashraf | Genomic prediction using phenotypes from pedigree lines with no markers[END_REF] were the first to investigate the impact on prediction accuracy when some wheat lines were not genotyped and only pedigree and phenotype information was available; the authors concluded that the ssGBLUP method for

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

4

deriving matrix H can provide higher prediction accuracy than either genomic or pedigree-based prediction. In plants, the ssGBLUP approach proposed by [START_REF] Ashraf | Genomic prediction using phenotypes from pedigree lines with no markers[END_REF] has been used

with a limited number of lines. The approach has not been tested on large datasets, e.g., every cropping cycle, CIMMYT's Global Wheat Program (GWP) generates thousands of new breeding lines that are candidates for field evaluation. Applying GS in the GWP is economically feasible

(1) when advancing breeding lines in the first preliminary yield trials to predict the performance of the selected lines in multi-environment trials or (2) for predicting a selected set of lines in different international target environments using as a training set the parents evaluated in Mexico and the progeny to be predicted in international environments such as those in South Asia.

In recent years, the GWP aimed at forming a large reference data set comprising 58,798 breeding lines, including their phenotypic and pedigree data from the last seven cropping cycles in Cd. Obregon (Mexico) and South Asia. This large reference set contains complete phenotypic data and pedigree information; however, only 29,484 of the lines have been genotyped.

Therefore, an H matrix that combines wheat lines having only molecular markers with those having pedigree and phenotype must be generated.

The main objectives of this study were (1) to use the large reference set for predicting the performance of wheat lines in several environments in South Asia; and (2) to perform the predictions using phenotypic, pedigree and genomic information to genetically evaluate the wheat lines using a single-step model that combines pedigree and marker information into a unified H matrix. Here we used information for genotyped and non-genotyped individuals combined, by applying the method proposed by [START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF] and [START_REF] Aguilar | Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF].

Prediction accuracy was studied using a G×E interaction multiplicative model (the reaction norm model of Jarquín et al., 2014) with pedigree information (A), genomic information (G), or both (H) and comparing its prediction accuracy results to those of a genomic model that does not include G×E interaction. This reaction norm model uses highly random dimensional matrices for the genomic and pedigree matrices. We also describe the statistical and computational challenges encountered when developing the pedigree and single-step models in such large data sets.

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

5

MATERIALS AND METHODS

Experimental data

The data set included a total of 58,798 wheat lines that were evaluated at the Norman E. (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016). Some of the lines were also evaluated under the same conditions in South Asia (Jalbapur, Ludhiana and Pusa in India; Faisalabad in Pakistan; and Jamalpur in Bangladesh)

during 2013-2016. Original data from each year comprise a large number of trials, each established using an alpha-lattice design with three replicates. The field management conditions under which each trial was established in each year are described in Table 1. The conditionslocation combinations will be referred to as environments. Table 2 shows the number of lines evaluated in each environment.

The basic model fitted to each of the 12 environments described in Table 2 comprises the random effects of the trials, the random effects of the replicates within trials, the random effects of the incomplete blocks within trials and replicates, and the random effects of the breeding lines.

A pedigree relationship matrix (A) for the 58,798 individuals was computed using a modified version of the software 'pedigreemm' [START_REF] Bates | pedigreemm: Pedigree-based mixed-effects models[END_REF]) that accounts for selfpollination; the latest version of the routines can be found at https://github.com/Rpedigree/pedigreeR. Given the dimensions of A, it is difficult to hold it in RAM memory and compute it. Appendix A shows the small R script (R Core Team, 2016) that was used to obtain and store the relationship matrix. It uses results from partitioned matrices to obtain the result and speed up the computations; R was recompiled from source and linked against OpenBLAS (http://www.openblas.net ). For further details on the computations, see Appendix A. In total, 29,484 individuals were genotyped using Genotyping by Sequencing (GBS) (e.g., [START_REF] Elshire | A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species[END_REF]. We kept all the SNP markers and imputed the missing values using observed data. Markers with minor allele frequency (MAF) of less than 0.05 were removed; after this process, 9,045 markers were available for prediction.

Statistical models

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

6

Recently, Jarquín et al. (2014) and [START_REF] López-Cruz | Increased prediction accuracy in wheat breeding trials using a marker × environment interaction genomic selection model[END_REF] proposed statistical models for performing genomic predictions taking into account G×E. The models were originally developed to incorporate genetic information from molecular markers and, in the case of Jarquín's model, it is also possible to incorporate environmental covariates. The Jarquín model has also shown to be useful when the genetic information is obtained from a pedigree (Pérez-Rodríguez et al., 2015). Here we describe Jarquín's model based on genomic and pedigree information. To speed up the computations and make them feasible, we re-parametrized the original model by using very well-known results from Cholesky decomposition and mixed models (e.g., [START_REF] Henderson | A simple method for computing the inverse of the numerator relationship matrix used in prediction of breeding values[END_REF][START_REF] Harville | Computational aspects of likelihood-based inference for variance components[END_REF].

Model 1: GE interaction using pedigree

The parametric GE interaction model takes into account the main effect of 𝐸 environments, the main effect of genotypes and the interaction between genotypes and the environment. In matrix notation, the model can be written as:

𝒚 = 𝜇𝟏 + 𝒁 𝐸 𝜷 𝐸 + 𝒁 𝑔 𝒖 1 + 𝒖 2 + 𝒆, (1)
where 𝒚 = (𝒚 1 , … , 𝒚 𝐸 )′ is the response vector, and 𝒚 𝑗 represents the observations in the j-th environment (𝑗 = 1, … , 𝐸). The general mean is 𝜇; 𝒁 𝐸 is an incidence matrix for environments, which is assumed multivariate with 𝜷 𝐸 ~𝑀𝑁(𝟎, 𝜎 𝐸 2 𝑰); 𝒁 𝑔 is an incidence matrix that connects genotypes with phenotypes; 𝒖 1 represents the random effect of genotypes; it is assumed multivariate with 𝒖 1 ~𝑀𝑁(𝟎, 𝜎 𝑢 2 𝑨); and 𝒖 2 represents the effect of GE interaction. We assume 𝒖 2 ~𝑀𝑁(𝟎, 𝜎 𝑔𝑒 2 (𝒁 𝑔 𝑮𝒁 𝑔 ′ )#(𝒁 𝐸 𝒁 𝐸 ′ )), where # denotes the Hadamard product (cell-by-cell) of the two matrices in parentheses (see Jarquín et al., 2014;Pérez-Rodríguez et al., 2015). Finally, we assume that the residuals are distributed as follows: 𝒆~𝑀𝑁(𝟎, 𝜎 𝑒 2 𝑰).

Since 𝑨 is positive definite and symmetric, it can be factored as 𝑨 = 𝑳𝑳′ by using Cholesky decomposition. Therefore, from (1):

𝒁 𝑔 𝒖 1 = 𝑑 𝒁 𝑔 𝑳𝒖 1 * (2)
where 𝒖 1 * ~𝑀𝑁(𝟎, 𝜎 𝑢 2 𝑰). Furthermore, it is not necessary to perform the 𝒁 𝑔 𝑳 product because for each row of the resulting matrix, we just need to copy the k-th row of L, where k is the column in the i-th row of 𝒁 𝑔 that is different from zero, that is, 𝒁 𝑔 (i, k) is equal to one. The matrix 𝒁 𝐸 𝒁 𝐸 ′ is block diagonal; blocks different from zero correspond to matrices with ones, i.e.,

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

7

𝒁 𝐸 𝒁 𝐸 ′ = ( 𝑱 1 𝑱 2 ⋱ 𝑱 𝐸 ), (3)
where 𝑱 𝑗 (𝑗 = 1, … , 𝐸) is a square matrix with ones whose dimensions correspond to the number of genotypes evaluated in environment 𝑗. Since 𝒁 𝐸 𝒁 𝐸 ′ is block diagonal, in order to compute 𝒁 𝑔 𝑨𝒁 𝑔 ′ #𝒁 𝐸 𝒁 𝐸 ′ , we just need to compute the corresponding block elements in the diagonal of

𝒁 𝑔 𝑨𝒁 𝑔 ′ = 𝒁 𝑔 𝑳𝑳 ′ 𝒁 𝑔 ′ . Let 𝒁 𝑔 𝑳 = 𝒁 ̃; then 𝑨 ̃= 𝒁 ̃𝒁 ̃′ = ( 𝒁 ̃11 𝒁 ̃12 ⋯ 𝒁 ̃1𝐸 𝒁 ̃21 𝒁 ̃22 ⋯ 𝒁 ̃2𝐸 ⋮ ⋮ ⋱ ⋮ 𝒁 ̃𝐸1 𝒁 ̃𝐸2 ⋯ 𝒁 ̃𝐸𝐸 ) ( 𝒁 ̃11 𝒁 ̃12 ⋯ 𝒁 ̃1𝐸 𝒁 ̃21 𝒁 ̃22 ⋯ 𝒁 ̃2𝐸 ⋮ ⋮ ⋱ ⋮ 𝒁 ̃𝐸1 𝒁 ̃𝐸2 ⋯ 𝒁 ̃𝐸𝐸 ) ′ ,
The block diagonal elements of 𝑨 ̃, can be computed as follows:

𝑨 ̃11 = ∑ 𝒁 ̃1𝑗 𝒁 ̃1𝑗 ′ 𝐸𝑛𝑣𝑖𝑟𝑜𝑛𝑚𝑒𝑛𝑡𝑠 = 𝑨 11 ⋮ (4) 𝑨 ̃𝐸𝐸 = ∑ 𝒁 ̃𝐸𝑗 𝒁 ̃𝐸𝑗 ′ 𝐸𝑛𝑣𝑖𝑟𝑜𝑛𝑚𝑒𝑛𝑡𝑠 = 𝑨 𝐸𝐸
where 𝑨 𝑗𝑗 , corresponds to the relationship matrix for individuals evaluated in environment 𝑗.

From (3) and (4) and by using Cholesky decomposition, the term 𝒁 𝑔 𝑨𝒁 𝑔 ′ #𝒁 𝐸 𝒁 𝐸 ′ can be obtained as follows:

𝒁 𝑔 𝑨𝒁 𝑔 ′ #𝒁 𝐸 𝒁 𝐸 ′ = 𝐵𝐷𝑖𝑎𝑔(𝑨 11 , … , 𝑨 𝐸𝐸 ) = 𝐵𝐷𝑖𝑎𝑔(𝑳 1 𝑳 1 ′ , … , 𝑳 𝐸 𝑳 𝐸 ′ ) = 𝑳 𝑔𝑒 𝑳 𝑔𝑒 ′ , ( 5
)
where 𝑳 𝑔𝑒 = 𝐵𝐷𝑖𝑎𝑔(𝑳 1 , … , 𝑳 𝐸 ). Therefore, from (5):

𝒖 2 = 𝑑 𝑳 𝑔𝑒 𝒖 2 * , ( 6
)
where 𝒖 2 * ~𝑀𝑁(𝟎, 𝜎 𝑔𝑒 2 𝑰) and "= 𝑑 " stands for equality in distribution.

Therefore, using results from (2) and (6), model (1) can be written as:

𝒚 = 𝜇𝟏 + 𝒁 𝐸 𝜷 𝐸 + 𝒁 𝑔 𝑳𝒖 1 * + 𝑳 𝑔𝑒 𝒖 2 * + 𝒆 (7)
Models ( 1) and ( 7) are equivalent, but model (7) has at least two advantages over model (1): (i) it avoids many matrix products, and (ii) it can be implemented relatively easily using the well known Gibbs sampler [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF] in the Bayesian framework.

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

8

Model 2: GE interaction using molecular markers

Let 𝑾 be a 𝑔 × 𝑝 matrix of standardized markers, where 𝑔 is the number of genotyped individuals and 𝑝 is the number of markers; let 𝑮 = 𝑾𝑾 ′ /𝑝 be the genomic relationship matrix [START_REF] López-Cruz | Increased prediction accuracy in wheat breeding trials using a marker × environment interaction genomic selection model[END_REF]. A model similar to (1) can be obtained by replacing 𝑨 with 𝑮.

Model 3: GE interaction using molecular markers and pedigree (single-step approach)

In this model, the information for genotyped and non-genotyped individuals is combined using the approach proposed by [START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF] and [START_REF] Aguilar | Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF]. A relationship matrix that includes full pedigree and genomic information is given by:

𝑯 = [ 𝑨 𝑛𝑛 + 𝑨 𝑔𝑛 ′ 𝑨 𝑔𝑔 -1 (𝑮 𝑎 -𝑨 𝑔𝑔 )𝑨 𝑔𝑔 -1 𝑨 𝑔𝑛 𝑨 𝑔𝑛 ′ 𝑨 𝑔𝑔 -1 𝑮 𝑎 𝑮 𝑎 𝑨 𝑔𝑔 -1 𝑨 𝑔𝑛 𝑮 𝑎 ],
where the matrix is divided based on whether the individuals have been genotyped or not.

Submatrices 𝑨 𝑔𝑔 , 𝑨 𝑛𝑛 and 𝑨 𝑔𝑛 are submatrices of 𝑨 containing relationships among genotyped individuals, among non-genotyped individuals and between genotyped and non-genotyped individuals, respectively [START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF][START_REF] Christensen | Single-step methods for genomic evaluation in pigs[END_REF]. 𝑮 𝑎 is an adjusted relationship matrix obtained from the genomic relationship matrix given by López-Cruz et al.

(2015), i.e., 𝑮 = 𝑾𝑾 ′ /𝑝 and 𝑨 𝑔𝑔 , that is:

𝑮 𝑎 = 𝛽𝑮 + 𝛼,
where 𝛽 and 𝛼 are obtained by solving the following system of equations:

𝐴𝑣𝑔(𝑑𝑖𝑎𝑔(𝑮))𝛽 + 𝛼 = 𝐴𝑣𝑔 (𝑑𝑖𝑎𝑔(𝑨 𝑔𝑔 )), 𝐴𝑣𝑔(𝑮)𝛽 + 𝛼 = 𝐴𝑣𝑔(𝑨 𝑔𝑔 ),
where 𝑮 𝑎 is a rescaled matrix such that: (i) the average of the diagonal elements is equal to the average of the diagonal elements of 𝑨 𝑔𝑔 , and (ii) the average of all the elements is equal to the average elements of 𝑨 𝑔𝑔 . See [START_REF] Christensen | Single-step methods for genomic evaluation in pigs[END_REF] for further details. Note that in this formulation based on H (and not its inverse), H does not need to be full rank.

Appendix B shows the R code that allows us to build the H matrix. A parametric GE interaction model takes into account the effect of 𝐸 environments, the main effect of genotypes and the interaction between genotypes and the environment. A model that uses information

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089 9 obtained from markers and pedigree can be obtained by replacing the A matrix in model ( 1) with the H matrix described above.

Model without GE interaction

Note that models that do not include the GE term can be derived from models 1-3 just by removing the corresponding random GE term. For example, by removing the term 𝒖 2

representing the effect of GE from model (1), it becomes

𝒚 = 𝜇𝟏 + 𝒁 𝐸 𝜷 𝐸 + 𝒁 𝑔 𝒖 1 + 𝒆
In this case, the resulting models are equivalent to the across-environment GBLUP model of [START_REF] López-Cruz | Increased prediction accuracy in wheat breeding trials using a marker × environment interaction genomic selection model[END_REF]. We include models without the GE term in order to compare the prediction accuracy of models with and without GE interaction. The single-environment model was not included because all the wheat lines included in the prediction of South Asian environments had complete pedigree and markers; thus developing matrix H for the single-step model does not make sense.

Assessing the model's predictive ability

The main interest of breeders is to predict the performance of non-evaluated lines in South Asian sites (Jalbapur, Ludhiana, Pusa in India; Faisalabad in Pakistan; and Jamalpur in Bangladesh). To mimic that situation, we designed a cross-validation scheme where we fitted the GE models (models 1-3) as well as models without GE using as the training set all available records under Drought, Late Heat, Optimal and Severe Drought conditions obtained in Cd. Obregon (Mexico), and 20% of available records in each of the South Asian sites assigned at random. In the prediction process, 80% of lines in the corresponding sites in the South Asian countries (India, Pakistan, and Bangladesh) were predicted using the rest of the records. A total of 20 random partitions (such as the ones described above) were generated.

The models' predictive ability was compared by using Pearson's correlation coefficient.

The models that use the A and H matrices included the phenotypic information of the 58,798 wheat lines, whereas the model that is based on markers only included information for 29,484 wheat lines that correspond to the individuals that were genotyped. The genotyped individuals are a subset of the individuals with pedigree information; therefore, lines in the testing set have pedigree and marker information. The numbers of individuals in the testing sets in South Asian

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089 10 sites are shown in Table 3, so in each random partition, the same individuals are predicted with three different models based on the A, G and H matrices.

Software

The models described above were fitted using a modified version of the BGLR package (de los Campos and Pérez-Rodríguez, 2015). The package was modified in order to accept as inputs big.matrix objects created using the bigmemory package [START_REF] Kane | Scalable Strategies for Computing with Massive Data[END_REF]. The bigmemory package was used to handle the huge matrices that have to be used during the analysis and also to take advantage of what in computer science is known as shared memory.

Once loaded into RAM memory, the data can be accessed from several processors, making it possible to perform a cross-validation relatively easily.

Data availability

The complete phenotypic and marker data can be found at http://genomics.cimmyt.org/wheat_50k/PG/.

RESULTS

Descriptive statistics

Figure 1 shows a boxplot of grain yield by location and median yield per location. From the plot it can be seen that the Optimal conditions had the highest grain yield, while the Late Heat and Severe Drought conditions had the worst grain yield. Yields in South Asian environments, especially in Pakistan and Bangladesh, were usually lower than in Mexican environments. Table 4 shows the number of lines evaluated in each environment, and lines in common between pairs of environments. It also shows sample correlations for grain yield for each pair of environments.

The number of lines evaluated in common between pairs of environments ranged from 537 to 4,735. The phenotypic sample correlation ranged from -0.05 to 0.53, which suggests large GE. Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089 11 entries of the H matrix are around 1.5, which stems from the standardization of G to be on the same scale as A.

Prediction accuracy of models without G×E

Table 4 shows the phenotypic correlations between pairs of environments. For example, the phenotypic correlation of the 4062 common wheat lines between environments B5I_OBR and B2I_OBR is 0.156, whereas the phenotypic correlation of the 1537 common wheat lines in B5I_OBR and STN_PUS (Pusa, India) is 0.210. In general, phenotypic correlations were not high, ranging from -0.051 to 0.481.

Table 5 shows the average Pearson's correlations between observed and predicted phenotypes and their corresponding standard deviation for the model without including G×E. The average correlations come from 20 random partitions with all the data records available in Mexico and 20% of the data available in South Asia. Note that these are the predictions of 80% of the entries included in the six South Asia environments. Prediction accuracies are relatively low, with those based on pedigree being slightly higher than those based on markers or on both pedigree and markers.

Prediction accuracy of G×E models

Table 6 shows the average of Pearson's correlations between observed and predicted phenotypes and its corresponding standard obtained using the same cross-validation scheme described above, but now including the G×E term. The predictive ability of models based on Pedigree, Markers and Pedigree + Markers is about the same, with pedigree prediction accuracy being higher than genomic and pedigree + genomic prediction accuracy in four environments (DLP_FAS, STN_JAM, STN_JBL, and STN_PUS). Ludhiana and Faisalabad under standard management conditions (0.3785, 0.2455, respectively) were the best predictive models for the genomic and pedigree + genomic model, respectively. Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

12 displays the correlation based on markers (G) versus that obtained based on pedigree (A),

where the prediction based on pedigree seems slightly better than that based on H (Fig. 3c).

Table 7 shows the percentage change in prediction accuracy of models with and without G×E. The % change was calculated as (𝑟 𝐺×𝐸 -𝑟 𝑛𝑜 𝐺×𝐸 )/𝑟 𝑛𝑜 𝐺×𝐸 × 100, where 𝑟 𝐺×𝐸 is the Pearson's correlation for a model with the G×E term and 𝑟 𝑛𝑜 𝐺×𝐸 is the Pearson's correlation for a model without the G×E term. From results in Table 7, it is clear that models that include the G×E term predict better than those that do not include G×E. For example, the G×E model using matrix H gave a 66% increase in prediction accuracy compared to the model using matrix H but without G×E. 8 presents the information needed to compute the weighted correlation for environment DLP_FAS; the columns present the information needed to compute the weighted correlation (note that this information was obtained from Table 4). The rest of the correlations were obtained using the approach described above. The gray bars represent the means of the correlations between observed and predicted values obtained from cross-validations. Note that in general the G×E models give good predictions, usually better than the phenotypic correlations. Although we are predicting 80% of the records in each of the South Asian environments, the correlations are higher than the phenotypical correlations between a given environment and the rest of the environments.

DISCUSSION

In wheat breeding, the cost of genotyping thousands of plants in segregating populations or in advanced generations makes the application of GS unfeasible. One possibility for solving this

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089 for assessing G×E; they found that the prediction accuracy of the single-step H matrix was higher than the accuracies achieved using the A and G matrices. The present study used seven selection cycles of CIMMYT wheat breeding with a total of 58,798 wheat lines evaluated in Cd.

Obregon and predicts several wheat lines in South Asian environments (India, Pakistan, and Bangladesh).

Genomic prediction accuracy for models with and without G×E

From the results in Tables 5, 6 and 7, it is clear that models that include the G×E term predict the environments in South Asia better than models that do not include the G×E term. The gain in

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089 14 prediction accuracy of models that include G×E ranges from 16 to 90% with an average of 40%.

However, models that do not incorporate G×E but use pedigree or high density molecular markers, or both, are still superior in terms of prediction accuracy than those that use only phenotypic data.

Genomic prediction accuracy versus phenotypic prediction accuracy of G×E models

In this study, we assessed the prediction accuracy of a large number of wheat lines evaluated in several environments and years in Cd. Obregon, Mexico, and predicted lines in several South Asian environments. For Ludhiana, Pusa, and Jabalpur, about 1227 wheat lines were predicted based on the performance of these lines in six environments in Cd. Obregon, plus the performance of about 57,000 wheat lines related to those to be predicted (1227) and evaluated in previous years in Cd. Obregon, Mexico.

Prediction accuracy was the correlation between the predicted values of the lines in Cd.

Obregon plus a low proportion of them (20%) in six environments in South Asia using three G×E models (A, G, and H) with the observed values of 80% of the lines in the six environments in South Asia (that were not phenotyped). The correlations for all the environments were around 0.25-0.27, except for Ludhiana in India, which showed higher prediction accuracy (0.36-0.37).

These genomic prediction accuracies were higher than the prediction accuracies computed from the common phenotypic correlations between all pairs of environments. These results indicated that the prediction accuracy with which breeders make selections in Cd. Obregon, Mexico, is lower than the accuracy they could obtain by performing genomic selection and prediction.

Although wheat breeders expect that lines selected in Cd. Obregon will perform well in South Asian environments, the results of this study should prompt them to increase research on genomic selection in Cd. Obregon (a very stable site with high radiation) of candidates for selection that will perform well in several environments in different South Asian countries (India, Pakistan, and Bangladesh).

The prediction accuracy of models with A, G, and H for models with or without G×E did not change much. This is an important result that allows, through the use of matrix H, using all phenotypic data to predict the genetic values of the unobserved wheat lines, thereby avoiding having to use only a subset of the phenotypes of those lines with pedigree and another subset of

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

15 phenotype data of lines with only marker data. Also the single-step method for computing H allows the inclusion of both components of the breeding value to be predicted, the parental average or between-family variability captured by the pedigree and Mendelian sample component (or with family variability) accounted by markers.

Big data used to derive pedigree and combine it with markers into the single-step prediction method with a G×E model

So far, no studies using plant breeding data on more than 50,000 lines have been reported in the genomic selection literature. This is the first study that shows that large training populations can provide genomic predictions that are more precise than phenotypic predictions. This is the first time that the theory used to develop and implement the pedigree system for such a large number of lines is reported in plant breeding. Although the models used for prediction are now well known, from the computational and statistical points of view, it is necessary to develop algorithms and data structures that allow researchers to handle the data and fit the models efficiently.

In this study, we used the G×E reaction norm model on a large data set in conjunction with pedigree, markers, or both, in genomic selection and prediction. We compared models including and excluding G×E and also in the genomic prediction literature there are plenty of examples where including those interactions significantly improved the prediction accuracy of untested individuals. The single-step method that combines the use of pedigree and markers through matrix H allows using all the available information. Also, the reaction norm G×E model allows borrowing information among positively correlated environments, although the predictive power of the model was similar to that of the model that includes markers only. Ashraf et al.

(2016) used the single-step H approach on a set of only 11,131 non-genotyped and 1176 genotyped wheat lines.

Animal breeders make extensive use of the fact that relationship matrix A has a very sparse inverse that can be computed directly from the pedigree, if all individuals (including those with no phenotype) are included [START_REF] Henderson | A simple method for computing the inverse of the numerator relationship matrix used in prediction of breeding values[END_REF][START_REF] Henderson | Best linear unbiased prediction of breeding values not in the model for records[END_REF]. This results in a sparse 𝑯 -1 structure

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089 16 as well [START_REF] Aguilar | Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF][START_REF] Christensen | Genomic prediction when some animals are not genotyped[END_REF], with a storage cost quadratic in the number of genotyped individuals and only linear in the number of non-genotyped individuals.

These sparse inverses exist for any level of autopolyploid species [START_REF] Kerr | Use of the numerator relationship matrix in genetic analysis of autopolyploid species[END_REF] and could potentially be used for prediction with large data sets. However, this would preclude the use of Cholesky decomposition used in (7).

CONCLUSIONS

This study shows how to solve statistical and computational challenges when incorporating and combining highly dimensional pedigree and genomic matrices into a single-step model for predicting unobserved individuals in other environments. We found that genomic prediction of genotyped and non-genotyped wheat lines produces higher prediction accuracy than that of lines predicted based on phenotypic data. The results provide evidence that the single-step approach that combines pedigree and marker information is useful for reducing genotyping costs while maintaining the prediction accuracy of unobserved individuals at relatively intermediate levels.

The incorporation of G×E models using a combination of pedigree and genomic information is another way of increasing the prediction accuracy of unobserved candidates for selection and offers plant breeders an important alternative for predicting germplasm evaluated under different environmental conditions.

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

21 Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

22 Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

23 Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

24 Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

25 Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

26 Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089 27 Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

29

APPENDIX

R script to obtain and store relationship matrix A

This script computes relationship matrix A.

Inputs:

1) A text file with pedigree information for the individuals that we are interested in. The file should have 3 columns separated by tabs, ID (Id of the individual), Sire and Dam.

2) A text file with the individuals that we are interested in.

Output: The relationship matrix.

To speed up the computations, we used dense partitioned matrixes and linked R against OpenBLAS (http://www.openblas.net). At the end of the process, the relationship matrix was also stored as a partitioned matrix on the hard disk in binary R format (RData). Below we detail the steps used to build the matrix.

Step 1: Read the data and compute the relationship matrix from the pedigree information cat("nrow=",nrow(a),"\n") cat("selfing=",sum(a[,1]==a[,2]),"\n") #Read the ids of individuals with phenotypic records ids=scan("GIDsForUSAIDprediction_20160406.csv")

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089 30 ids=as.character(ids) pede=editPed(sire=a$SIRE,dam=a$DAM,label=a$ID,verbose=TRUE) ped=with(pede, pedigree(label=label, sire=sire, dam=dam))

Now use the relfactor function for the pedigree, that is

𝑨 𝒇𝒖𝒍𝒍 = 𝑿 𝒇𝒖𝒍𝒍 ′𝑿 𝒇𝒖𝒍𝒍
where 𝑿 𝒇𝒖𝒍𝒍 is an upper triangular, sparse (right) Cholesky factor of the relationship matrix. In this case, 𝑿 𝒇𝒖𝒍𝒍 is a matrix with n=177,376 rows and the same number of columns. The code for obtaining the relfactor is given below.

Xfull=relfactor(ped)

We do not need 𝑨 𝒇𝒖𝒍𝒍 ; we just need a subset of this matrix with the 58,798 individuals so we can take a subset of 58,798 columns from 𝑿 𝒇𝒖𝒍𝒍 . The columns correspond to the individuals that we are interested in. Let X be the resulting matrix; then

𝑨 = 𝑿′𝑿

where X has n=177,376 rows and p=58,798 columns. The R code for obtaining this matrix is shown below.

index=ped@label%in%ids X=Xfull[,index]

Step 2: Partition the relationship factor

Since X is a huge matrix, it is very difficult to obtain 𝑨 directly; furthermore, since X is sparse, the product cannot be parallelized easily. So we partitioned X into several sub-matrices and saved the sub-matrices as binary files that can later be retrieved in order to perform the product.

For example:

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2 where 𝑿 𝑖𝑗 is a sub-matrix obtained from 𝑿.

The R code below was used to partition matrix X into 5 sub-matrices and save the results to binary files. cat("Submatrix: ",i," ",j,"\n"); cat("from_row: ",from_row,"\n");

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089 32 cat("to_row: ",to_row,"\n"); cat("from_column: ",from_column,"\n");

#Conventional matrix object so that we can use #optimized dense matrix products Xij=as.matrix(X[from_row:to_row,from_column:to_column])

save(Xij,file=paste("X_",i,j,".RData",sep="")) } }

Step 3: Compute the relationship matrix using the partitioned matrices from step 2

Given the partition of the relationship factor, we can compute the A matrix as follows: # a problem here is to divide A neatly between genotyped and not genotyped individuals. # Usually we use sparse operators and this is easier. # here I use the colnames and should be efficient # in principle H is (SEMI-)positive definite but can be quite bad conditioned, # e.g. if there are pedigree errors or label switching save(H,file="H.Rdata")

𝑨 = 𝑿′𝑿 = (
# ---------------------------------------# # option 1 to construct H not its inverse # ---------------------------------------# # use
Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

36 Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

37 Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

39 Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

40 6 using the H matrix. Black bars represent a weighted mean of the phenotypic correlation of a given environment and the rest of environments in Table 4; for example, for DPL_FAS, the weighted correlation can be obtained using data shown in Table 8. FAS=Faisalabad, Pakistán; JAM=Jamalpur, Bangladesh; JBL=Jabalpur, India; LDH=Ludhiana, India, OBR=Obregon, Mexico; PUS=Pusa, India. Standard management conditions (STN), and delayed planting conditions (DLP).

Figure 2

 2 Figure 2 displays the distribution of the diagonal entries for matrices A, H and G. Note that in the A matrix, the diagonal entries are around ~1.5; in this case, 𝑎(𝑖, 𝑖) = 1 + 𝐹 𝑖 , where 𝐹 𝑖 is the inbreeding coefficient of the i-th individual. The diagonal entries of the G matrix are around 1.0, reflecting the fact that the markers were centered and standardized. The diagonal

Figure 3

 3 Figure 3(a)-(c) shows scatterplots of predictive correlations for each of 20 crossvalidations across the six environments in South Asia. Figure 3a depicts the correlations between predicted values based on markers (G) versus those based on matrix H and shows that prediction accuracy based on G was superior to that obtained based on H. Figure 3b

Figure 4

 4 Figure4presents a bar plot of correlations for each predicted environment in South Asia using the H matrix. Bars in black represent the mean of the weighted phenotypic correlation of a given environment and the rest of the environments in Table4. The phenotypic correlation for environment j in South Asia can be obtained as follows: 𝑟 𝑗 = ∑

13

  problem would be to augment the numerical relationship matrix (A) of all individuals with the genomic relationship matrix (G) of the genotyped individuals and perform predictions based on the resulting complete H, which will allow performing prediction of non-genotyped individuals in the testing set. Augmenting matrix A by using only a fraction of the genotyped individuals would reduce genotyping costs. Furthermore, as described by[START_REF] Christensen | Single-step methods for genomic evaluation in pigs[END_REF], the single-step method allows the genomic relationship matrix of genotyped individuals to be extended using pedigree information to a combined relationship matrix H of all individual plants or lines. This allows using all phenotypic data and not only those phenotypes that have pedigree and marker information; this extra phenotypic information should also enhance prediction accuracy. This makes the models and methods developed by[START_REF] Misztal | Computing procedures for genetic evaluation including phenotypic, full pedigree, and genomic information[END_REF],[START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF] and[START_REF] Aguilar | Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF] 2011) very attractive for predicting unobserved and nongenotyped plants.In a recent article,[START_REF] Fernando | A class of Bayesian methods to combine large numbers of genotyped and non-genotyped animals for whole genome analyses[END_REF] proposed a single-step Bayesian regression strategy that allows using all genotyped and non-genotyped individuals by means of imputed marker covariates for non-genotyped individuals. The advantage of the Bayesian approach over the single-step BLUP is that it does not require computing the inverse of G. However, this modelhas not yet been applied to realistic datasets.The single-step approach of[START_REF] Misztal | Computing procedures for genetic evaluation including phenotypic, full pedigree, and genomic information[END_REF],[START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF] and[START_REF] Aguilar | Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF] 2011) was used in dairy cattle studies and first applied to plant breeding data by[START_REF] Ashraf | Genomic prediction using phenotypes from pedigree lines with no markers[END_REF] in a set of 1, 176 genotyped CIMMYT wheat lines and 11,131 non-genotyped wheat lines tested in five environments in Cd. Obregon, Mexico, during the 2012-2013 cycle. The authors developed optimized weighting factors for matrix H and applied a multivariate method

  planting and 5 irrigations (B5I); Bed planting and 2 irrigations (B2I); Zero till, bed planting and 5 irrigations (Z5I); Zero till, bed planting and 2 irrigations (Z2I); Melgas flat planting (MEL); Drip irrigation to impose drought in flat (DRM); Drip irrigation to impose drought, beds (DRB); Early heat (EHT); Late heat (LHT).

  ***********************\n")

Figure 1 :

 1 Figure 1

Figure 2 :

 2 Figure 2

  Figure 3a

Figure 3 :

 3 Figure 3c

Figure 4 DLP_FASFigure 4 :

 44 Figure 4

Table 1 :

 1 Description of the conditions under which the 58,798 wheat lines were evaluated in different environments.

	Description	Field Management
		Conditions
	Standard management conditions (STN)	Optimal
	Delayed planting (DLP)	Late Heat
	Bed planting and 5 irrigations (B5I)	Optimal
	Bed planting and 2 irrigations (B2I)	Drought
	Zero till, bed planting and 5 irrigations (Z5I)	Optimal
	Zero till, bed planting and 2 irrigations (Z2I)	Drought
	Melgas flat planting and 5 irrigations (MEL)	Optimal
	Melgas flat planting and drip irrigation (DRM) Severe Drought
	Bed planting and drip irrigation (BDR)	Severe Drought
	Early heat (EHT)	Early Heat
	Late heat (LHT)	Late Heat

Table 2 :

 2 Number of lines evaluated in different environments during 2009-2016 by the GWP.

	Environment	Number of lines evaluated
	B5I_OBR	56,964
	B2I_OBR	4,063
	DRB_OBR	5,913
	EHT_OBR	2,188
	LHT_OBR	4,736
	MEL_OBR	4,735
	DLP_FAS	1,547
	STN_FAS	1,547
	STN_JAM	537
	STN_JBL	1,548
	STN_LDH	1,548
	STN_PUS	1,548
	FAS=Faisalabad, Pakistan; JAM=Jamalpur,
	Bangladesh; JBL=Jabalpur, India;
	LDH=Ludhiana, India, OBR=Cd, Obregon,
	Mexico; PUS=Pusa, India. Standard management
	conditions (STN); Delayed planting (DLP); Bed

Table 3 :

 3 Number of individuals in the testing set in South Asian sites.

	Environment	Number of individuals in the
		testing set
	DLP_FAS	1237
	STN_FAS	1237
	STN_JAM	429
	STN_JBL	1238
	STN_LDH	1238
	STN_PUS	1238

Table 4 :

 4 Number of genotypes (diagonal and upper triangular) and sample phenotypic correlations (lower triangular) by environment.

	Env	B5I_OBR B2I_OBR DRB_OBR EHT_OBR LHT_OBR MEL_OBR DLP_FAS STN_FAS STN_JAM STN_JBL STN_LDH STN_PUS
	B5I_OBR	56964	4062	4090	2187	4734	4735	1537	1537	532	1537	1537	1537
	B2I_OBR	0.156	4063	4063	2186	4063	4062	1515	1515	530	1515	1515	1515
	DRB_OBR -0.050	0.534	5913	2186	4091	4090	1535	1535	530	1535	1535	1535
	EHT_OBR	0.479	0.186	-0.051	2188	2187	2187	1062	1062	532	1062	1062	1062
	LHT_OBR	0.203	0.262	0.167	0.199	4736	4734	1537	1537	532	1537	1537	1537
	MEL_OBR	0.370	0.238	0.117	0.354	0.169	4735	1537	1537	532	1537	1537	1537
	DLP_FAS	0.154	0.094	0.111	0.131	0.067	0.174	1547	1547	537	1547	1547	1547
	STN_FAS	0.124	0.120	0.167	0.009	0.029	0.102	0.338	1547	537	1547	1547	1547
	STN_JAM	0.228	0.146	0.130	0.160	0.079	0.113	0.170	0.206	537	537	537	537
	STN_JBL	0.188	0.176	0.168	0.082	0.136	0.143	0.235	0.263	0.136	1548	1548	1548
	STN_LDH	0.225	0.079	0.078	0.190	0.040	0.168	0.206	0.286	0.382 0.140	1548	1548
	STN_PUS	0.210	0.137	0.099	0.117	0.025	0.173	0.280	0.241	0.481 0.255	0.222	1548
	FAS=Faisalabad, Pakistan; JAM=Jamalpur, Bangladesh; JBL=Jabalpur, India; LDH=Ludhiana,
	India, OBR=Obregon, Mexico; PUS=Pusa, India. Standard management conditions (STN);
	Delayed planting (DLP); Bed planting and 5 irrigations (B5I); Bed planting and 2 irrigations

(B2I); Zero till, bed planting and 5 irrigations (Z5I); Zero till, bed planting and 2 irrigations (Z2I); Melgas flat planting (MEL); Drip irrigation to impose drought in flat (DRM); Drip irrigation to impose drought, beds (DRB); Early heat (EHT); Late heat (LHT).

Table 5 :

 5 Correlations between predicted and observed values obtained using the cross-validation where all the wheat lines from Cd. Obregon, Mexico, plus 20% of the wheat lines in each of the environments in India, Pakistan, and Bangladesh were used in the training set to predict 80% of the lines in the corresponding environments in India, Pakistan, and Bangladesh. The highest correlations in each environment are in boldface.

			Model without G×E
	Environment	Pedigree (A)	Markers (G)	Pedigree + Markers (H)
	DLP_FAS	0.2113 (0.0304) 0.1716 (0.0104)	0.1834 (0.0135)
	STN_FAS	0.1611 (0.0181) 0.1235 (0.0129)	0.1455 (0.0120)
	STN_JAM	0.2448 (0.0251) 0.1861 (0.0189)	0.1992 (0.0213)
	STN_JBL	0.2480 (0.0184) 0.1928 (0.0154)	0.2075 (0.0163)
	STN_LDH	0.2554 (0.0158) 0.2472 (0.0104)	0.2477 (0.0094)
	STN_PUS	0.2361 (0.0143) 0.1989 (0.0112)	0.2117 (0.0107)
	FAS=Faisalabad, Pakistan; JAM=Jamalpur, Bangladesh; JBL=Jabalpur, India;
	LDH=Ludhiana, India, OBR=Obregon, Mexico; PUS=Pusa, India. Standard
	management conditions (STN), and delayed planting conditions (DLP).

Table 6 :

 6 Correlations between predicted and observed values obtained using the cross-validation where all the wheat lines from Cd. Obregon, Mexico, plus 20% of the wheat lines in sites in India, Pakistan, and Bangladesh were used in the training set to predict 80% of the lines in the corresponding sites in India, Pakistan, and Bangladesh. The highest correlations in each environment are in boldface.

			G×E model
	Environment	Pedigree (A)	Markers (G)	Pedigree + Markers (H)
	DLP_FAS	0.2462 (0.0294) 0.2327 (0.0132)	0.2345 (0.0123)
	STN_FAS	0.2360 (0.0227) 0.2414 (0.0180)	0.2455 (0.0175)
	STN_JAM	0.2942 (0.0414) 0.2681 (0.0293)	0.2656 (0.0309)
	STN_JBL	0.2921 (0.0183) 0.2741 (0.0163)	0.2739 (0.0165)
	STN_LDH	0.3699 (0.0109) 0.3785 (0.0157)	0.3651 (0.0155)
	STN_PUS	0.2842 (0.0175) 0.2622 (0.0191)	0.2684 (0.0185)
	FAS=Faisalabad, Pakistan; JAM=Jamalpur, Bangladesh; JBL=Jabalpur, India;
	LDH=Ludhiana, India, OBR=Obregon, Mexico; PUS=Pusa, India. Standard
	management conditions (STN), and delayed planting conditions (DLP).

Table 7 .

 7 Comparing the predictive ability of models with and without G×E. (𝑟 𝐺×𝐸 -𝑟 𝑛𝑜 𝐺×𝐸,𝐶 )/𝑟 𝑛𝑜 𝐺×𝐸 × 100, where 𝑟 𝐺×𝐸 is the Pearson's correlation for a model with the G×E term and 𝑟 𝑛𝑜 𝐺×𝐸 is the Pearson's correlation for a model without the G×E term.

	Environment	% Change Pedigree (A) Markers (G) Pedigree + Markers (H)
	DLP_FAS	16.52	35.61	26.88
	STN_FAS	46.49	95.47	65.91
	STN_JAM	20.18	44.06	34.59
	STN_JBL	17.78	42.17	32.10
	STN_LDH	44.83	53.11	52.81
	STN_PUS	20.37	31.83	23.85
	% Change = FAS=Faisalabad, Pakistan; JAM=Jamalpur, Bangladesh; JBL=Jabalpur, India;
	LDH=Ludhiana, India, OBR=Obregon, Mexico; PUS=Pusa, India. Standard
	management conditions (STN), and delayed planting conditions (DLP).

Table 8 :

 8 Phenotypic correlations and numbers of lines in common between DLP_FAS and the rest of the environments in Mexico and South Asia.

	j		k		𝑟 𝑗𝑘	𝑛 𝑗𝑘	𝑛 𝑗𝑘 /𝑛 𝑗 𝑛 𝑗𝑘 /𝑛 𝑗 𝑟 𝑗𝑘
		Env. in South Asia		Other environments
	1	DLP_FAS	1	B5I_OBR	0.154 1537 0.099	0.015
	1	DLP_FAS	2	B2I_OBR	0.094 1515 0.098	0.009
	1	DLP_FAS	3	DRB_OBR	0.111 1535 0.099	0.011
	1	DLP_FAS	4	EHT_OBR	0.131 1062 0.069	0.009
	1	DLP_FAS	5	LHT_OBR	0.067 1537 0.099	0.006
	1	DLP_FAS	6	MEL_OBR	0.174 1537 0.099	0.017
	1	DLP_FAS	7	STN_FAS	0.338 1547 0.100	0.033
	1	DLP_FAS	8	STN_JAM	0.170	537 0.035	0.005
	1	DLP_FAS	9	STN_JBL	0.235 1547 0.100	0.023
	1	DLP_FAS	10	STN_LDH	0.206 1547 0.100	0.020
	1	DLP_FAS	11	STN_PUS	0.280 1547 0.100	0.028
					𝑛 1 = 15448		𝑟 1 = 0.18

  ). , DOI : 10.3835/plantgenome2016.09.0089

		𝑿 11 𝑿 12
	𝑿 =	𝑿 21 𝑿 22 𝑿 31 𝑿 32 𝑿 41 𝑿 42	; 𝑿 ′ = ( 𝑿 11 ′ 𝑿 12 ′	𝑿 21 ′ 𝑿 22 ′	𝑿 31 ′ 𝑿 32 ′	𝑿 41 ′ 𝑿 42 ′	′ 𝑿 51 𝑿 52 ′ )
		( 𝑿 51 𝑿 52)
			31

  ,genotyped]=A[nongenotyped,genotyped]%*%(Aggi%*%G) #H[genotyped,nongenotyped]=G%*%Aggi%*%A[genotyped,nongenotyped] H[genotyped,nongenotyped]=t(H[nongenotyped,genotyped]) H[nongenotyped,nongenotyped]=A[nongenotyped,nongenotyped] + A[nongenotyped,genotyped]%*%(Aggi%*%(G-Agg)%*%Aggi)%*%A[genotyped,nongenotyped] cat(mean(diag(H)),mean(H),"\n")

	expression (4) in Legarra-Aguilar-Misztal 2009
	H=matrix(NA,ncol(A),nrow(A))
	colnames(H)=colnames(A)
	rownames(H)=rownames(A)
	H[genotyped,genotyped]=G
	H[nongenotyped

Plant Gen. Accepted Paper, posted 01/25/2017. doi:10.3835/plantgenome2016.09.0089

ACKNOWLEDGMENTS

We are grateful to field and lab workers of CIMMYT's Global Wheat Program and to national program researchers who collected the data used in this study. This work was funded by USAID Feed the Future Innovation Lab for applied Wheat Genomics with headquarters at Kansas State University in collaboration with CIMMYT and their partners in South Asia, specifically in India and Pakistan. We also thank the Bill and Melinda Gates Foundation for providing financial support for genotyping the materials used here.

Version postprint

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

Version postprint

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

Version postprint

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. ( 2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089

Version postprint

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089 33 for(i in 1:n_submatrix) { cat("i=",i,"\n") load(paste("X_",i,"1.RData",sep="")) A11=A11+crossprod(Xij); } save(A11,file="A11.RData")

The rest of the matrices can be similarly obtained. With this approach and using 8 cores for the matrix product, we obtained the 58, 798 × 58, 798 matrix A in less than 3 hours in the CIMMYT-BSU server (which has 12 Intel (R) Xeon Cores @ 3.47 GHz and ~ 48 Gb of RAM).

R script to obtain the H matrix

The script presented below computes a relationship matrix H including full pedigree and genomic information (see equation 4 in [START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF]. It adjusts the elements of genomic relationship matrix G, so that the entries of relationship matrix A share the same scale [START_REF] Christensen | Single-step methods for genomic evaluation in pigs[END_REF].

Inputs:

1) A and G matrices. The row and column names of both matrices include the IDs of the individuals.

Output:

1) H matrix.

#Clean workspace rm(list=ls())

#Load A load("../output/A11.RData") load("../output/A12.RData") load("../output/A21.RData") load("../output/A22.RData")

Version postprint

Comment citer ce document : Pérez-Rodríguez, P., Crossa, J., Rutkoski, J., Poland, J., Singh, R., Legarra, A., Autrique, E., De Los Campos, G., Burgueño, J., Dreisigacker, S. (2017). Single-step genomic and pedigree genotype x environment interaction models for predicting wheat lines in international environments. Plant Genome, 10 (2). , DOI : 10.3835/plantgenome2016.09.0089


Sophie Tronnet

Christophe Garcie

Alexander O Brachmann

Jörn Piel

Eric Oswald
email: eric.oswald@inserm.fr

Patricia Martin

High iron supply inhibits the synthesis of the genotoxin colibactin by pathogenic Escherichia coli through a non-canonical Fur/RyhB-mediated pathway

Keywords: Iron, colibactin, RyhB, Fur

The genotoxin colibactin is a secondary metabolite produced by a variety of pathogenic Enterobacteria, and is associated with colon cancer development and acute systemic infections. The colibactin biosynthesis requires the enzymatic activity of the phosphopantetheinyl transferase ClbA. We recently evidenced that two master regulators of bacterial iron homeostasis, i.e. the ferric uptake regulator (Fur) and the small regulatory noncoding RNA RyhB, were involved in the regulation of the clbA transcription and of the colibactin production. In this study we investigated the impact of high iron supply on clbA transcription and colibactin production in wild type, ΔryhB, Δfur and ΔryhB Δfur strains. This revealed that high iron resulted in decreased synthesis of the genotoxin colibactin through both pathways dependent and independent of Fur/RyhB. This work highlights the complex regulatory mechanism that controls an important bacterial virulence and carcinogenesis factor by regulators of bacterial iron homeostasis.

Introduction

The regulation of iron homeostasis in bacteria mainly relies on ferric uptake regulator (Fur) protein and the small regulatory non-coding RNA (sRNA) RyhB [START_REF] Massé | Small RNAs controlling iron metabolism[END_REF]Troxell and Hassan 2013;[START_REF] Porcheron | Interplay between iron homeostasis and virulence Fur and RyhB as major regulators of bacterial pathogenicity[END_REF]. Fur is a global transcriptional regulator that controls the transcription of over 90 genes involved in iron uptake, storage and metabolism. Fur is a well-known repressor. However, a growing literature reports Fur as an activator (Yu and Genco 2012;Troxell and Hassan 2013). Fur has been shown to play an important role in the modulation of bacterial virulence (Gancz et al. 2006;[START_REF] Payne | Iron and pathogenesis of Shigella iron acquisition in the intracellular environment[END_REF]Ledala et al. 2010;Torres et al. 2010;Tobe et al. 2014). Bacterial small regulatory RNA (sRNAs) are commonly known to repress gene expression by base pairing to target mRNAs. In contrast, a growing number of examples of translational activation and mRNA stabilization by sRNAs have now Comment citer ce document : Tronnet, S., Garcie, C., Brachmann, A. O., Piel, Oswald, E., Martin, P. (2017). High iron supply inhibits the synthesis of the genotoxin colibactin by pathogenic Escherichia coli through a non-canonical Fur/RyhB-mediated pathway. Pathogens and Disease, 75 (5). , DOI : 10.1093/femspd/ftx066 3 been documented [START_REF] Papenfort | Target activation by regulatory RNAs in bacteria[END_REF]. Although Fur and iron directly regulate sRNA RyhB [START_REF] Massé | Effect of RyhB small RNA on global iron use in Escherichia coli[END_REF][START_REF] Sobrero | The bacterial protein Hfq much more than a mere RNA-binding factor[END_REF], strong evidence suggests that RyhB itself is sufficient to influence the level of free intracellular iron [START_REF] Massé | Effect of RyhB small RNA on global iron use in Escherichia coli[END_REF].

RyhB was shown to be implicated in virulence-associated processes in pathogenic bacteria such as Shigella flexneri [START_REF] Oglesby | Fur regulates acid resistance in Shigella flexneri via RyhB and ydeP[END_REF], S. dysenteriae [START_REF] Murphy | RyhB an Iron-Responsive Small RNA Molecule Regulates Shigella dysenteriae Virulence[END_REF], Listeria monocytogenes [START_REF] Sesto | Bacterial and cellular RNAs at work during Listeria infection[END_REF] or Escherichia coli [START_REF] Porcheron | The small RNA RyhB contributes to siderophore production and virulence of uropathogenic Escherichia coli[END_REF].

It is observed over the last 30 years a change in the composition of commensal E. coli population in industrialized countries, with an increase in the phylogenetic group B2 (Tenaillon et al. 2010;[START_REF] Massot | Phylogenetic virulence and antibiotic resistance characteristics of commensal strain populations of Escherichia coli from community subjects in the Paris area in 2010 and evolution over 30 years[END_REF]. Genetic analysis of different E. coli strains showed that those of the phylogenetic group B2 have developed a greater ability to acquire iron [START_REF] Schubert | Role of intraspecies recombination in the spread of pathogenicity islands within the Escherichia coli species[END_REF]. Colibactin is a toxin synthesized by a diversity of pathogenic Enterobacteria, including approximately 50% of the E. coli stains that belong to the phylogenetic group B2 [START_REF] Nougayrède | Escherichia coli induces DNA double-strand breaks in eukaryotic cells[END_REF]. This genotoxin induces DNA double strand breaks, senescence, chromosomal abnormalities in enterocytes both in vitro and in vivo (Cuevas-Ramos et al. 2010;[START_REF] Secher | Escherichia coli producing colibactin triggers premature and transmissible senescence in mammalian cells[END_REF][START_REF] Payros | Maternally acquired genotoxic Escherichia coli alters offspring's intestinal homeostasis[END_REF], and constitutes a risk factor for the onset of colon tumor genesis [START_REF] Arthur | Intestinal inflammation targets cancerinducing activity of the microbiota[END_REF]. Colibactin also exacerbates the development of systemic pathologies such as sepsis [START_REF] Marcq | The genotoxin colibactin exacerbates lymphopenia and decreases survival rate in mice infected with septicemic Escherichia coli[END_REF]) and neonatal meningitis [START_REF] Mccarthy | The genotoxin Colibactin is a determinant of virulence in Escherichia coli K1 experimental neonatal systemic infection[END_REF]. The synthesis of colibactin relies on the pks gene cluster [START_REF] Nougayrède | Escherichia coli induces DNA double-strand breaks in eukaryotic cells[END_REF]) and requires the enzymatic activity of the phosphopantetheinyl transferase ClbA that is encoded on the pks island. Recently we showed that the transcription of clbA and the production of colibactin were regulated by Fur and RyhB (Tronnet et al. 2016). Therefore, iron could constitute a key environmental factor impacting the virulence of E. coli strains producing colibactin.

Here, we addressed the question of the role of high iron supply on the production of colibactin. We demonstrate that high iron inhibits the production of colibactin genotoxin Comment citer ce document : Tronnet, S., Garcie, C., Brachmann, A. O., Piel, Oswald, E., Martin, P. (2017). High iron supply inhibits the synthesis of the genotoxin colibactin by pathogenic Escherichia coli through a non-canonical Fur/RyhB-mediated pathway. Pathogens and Disease, 75 (5). , DOI : 10.1093/femspd/ftx066 4 through a complex regulatory mechanism that involves Fur, RhyB and an additional regulator of bacterial iron homeostasis.

Materials and Methods

Bacterial strains, mutagenesis and growth conditions

Bacterial strains used in this study are listed in Table 1. For genetic manipulations, E. coli strains were routinely grown at 37°C under shaking in 3ml of Lennox L broth (LB, Invitrogen). Kanamycin (50 μg/ml) or chloramphenicol (25 μg/ml) were added to the medium when required. Inactivation of E. coli genes was undertaken using phage λ Red recombinase (Datsenko and Wanner 2000). The primers used to construct mutants are shown in Table 1.

Allelic exchanges were confirmed by PCR. For megalocytosis assay, genotoxicity assay and bioluminescence measurements, E. coli strains were grown overnight in DMEM containing Hepes (DMEM-Hepes, Gibco, 0.034 µM Fe) supplemented with 100 µM FeCl 3 (20.66 µM Fe) when required, at 37°C with shaking. The overnight cultures were then diluted 1:50 in DMEM-Hepes supplemented or not with FeCl 3 and grown until OD 600nm = 0.6.

Luciferase measurements

Promoter activity of gene clbA was determined by time-course quantification of luciferase as previously described (Tronnet et al. 2016). E. coli was grown overnight in DMEM-Hepes, and subcultured into DMEM-Hepes or DMEM-Hepes supplemented with FeCl 3 up to OD 600nm = 0.6. Hundred microliter samples were then used to inoculate black 96-well plate (Greiner Bio-one) and grown without shaking at 37°C in a luminometer (Tecan Infinite Pro microplate reader). Both OD 600nm and light emission (relative light units, RLU, 6000ms aperture/sample) were recorded every 30 min, simultaneously.

Comment citer ce document : Tronnet, S., Garcie, C., Brachmann, A. O., Piel, Oswald, E., Martin, P. (2017). High iron supply inhibits the synthesis of the genotoxin colibactin by pathogenic Escherichia coli through a non-canonical Fur/RyhB-mediated pathway. Pathogens and Disease, 75 (5). , DOI : 10.1093/femspd/ftx066 5 N-myristoyl-D-asparagine quantification by liquid chromatography/mass spectrometry Quantification of colibactin prodrug motif was performed as previously described (Garcie et al. 2016). Briefly, E. coli strain M1/5 was grown in DMEM-Hepes with or without FeCl 3 at 37°C for 18h. Supernatants of cultures were obtained by centrifugation and filtered.

Quantification experiments were conducted with ultra-performance liquid chromatographyhigh resolution/heated electrospray ionization mass spectrometry (UPLC-HR/HESI-MS). The data were recorded on a Thermo Scientific Q Exactive Hybrid Quadrupole-Orbitrap Mass Spectrometer coupled to a Dionex Ultimate 3000 UPLC. The MS was operated in positive ionization mode at a scan range of 200-500 m/z and a resolution of 35,000. The spray voltage was set to 3.5 kV, the S-lens to 35, the auxiliary gas heater temperature to 438°C and the capillary temperature to 270°C. Absolute quantification was achieved by using a Schotten-Baumann reaction-derived N-myristoyl-L-asparagine (isomer of the N-myristoyl-D-asparagine colibactin cleavage product) as a standard. Data were obtained from undiluted cell free sample supernatants, concentrations were calculated using Thermo Xcalibur 2.2 Quan Browser.

Megalocytosis assay

Colibactin induces megalocytosis in cultured eukaryotic cells, manifest as progressive enlargement of the cell body and nucleus and the absence of mitosis. We quantified the extent of colibactin-induced megalocytosis using a methylene blue binding assay [START_REF] Martin | Interplay between siderophores and colibactin genotoxin biosynthetic pathways in Escherichia coli[END_REF]). E. coli strains were added to HeLa cells at multiplicities of infection (MOIs) of 200 and/or 100, co-cultured for 4h and washed. Cells were then incubated for 72h with cell culture medium containing 200 μg/mL gentamicin followed by staining with methylene blue. Dye binding was determined spectrophotometrically at OD 660 .

Comment citer ce document : Tronnet, S., Garcie, C., Brachmann, A. O., Piel, Oswald, E., Martin, P. (2017). High iron supply inhibits the synthesis of the genotoxin colibactin by pathogenic Escherichia coli through a non-canonical Fur/RyhB-mediated pathway. Pathogens and Disease, 75 (5). , DOI : 10.1093/femspd/ftx066 6

Genotoxicity assay

The capacity of colibactin to engender double strand DNA breaks was determined in HeLa cells by γ-H2AX immunofluorescence analysis [START_REF] Nougayrède | Escherichia coli induces DNA double-strand breaks in eukaryotic cells[END_REF]; this assay monitors the phosphorylation of histone H2AX, a sensitive marker of ds DNA breaks. HeLa cells (1.5 × 10 4 cells in 200 μl Dulbecco's modified Eagle's medium supplemented with 10% fetal calf serum, 80 µg/mL gentamicin and 0.1 unit/ml bovine insulin) were dispensed into 96-well cell culture plates. After incubation at 37 o C in a 5% CO 2 atmosphere for 24h, the cells were washed and incubated with bacteria at MOI = 50 bacteria per cell. After 4h infection at 37 o C in 5% CO 2 , cells were washed 3 times with Hank's Balanced Salts Solution and incubated at 37 o C in cell culture medium for 3h with 200 μg/ml gentamicin. Cells were fixed in the plate with 4% paraformaldehyde and processed as previously described [START_REF] Martin | Interplay between siderophores and colibactin genotoxin biosynthetic pathways in Escherichia coli[END_REF]. Rabbit monoclonal anti-γ-H2AX antibody #9718 (Cell Signaling Technology Inc, Danvers MA) was diluted 1:500 in blocking solution and incubated for 2h at room temperature. IRDyeTM 800CW-conjugated goat anti-rabbit secondary antibody was diluted 1:100 in blocking solution and incubated for 1 h. RedDot2 (Biotium) was used for DNA labelling. DNA and γ-H2AX were visualized using an Odyssey Infrared Imaging Scanner (Li-Cor ScienceTec, Les Ulis, France) using 680nm and 800nm channels for RedDot2 and IRDyeTM800 respectively.

Relative fluorescence units for γ-H2AX per cell (as determined by γ-H2AX divided by DNA content) were divided by vehicle controls to determine percentage change in phosphorylation of H2AX levels relative to control. All experiments were carried out in triplicate.

Statistical analysis

Statistical analyses were conducted using GraphPad Prism 6.0. The mean with standard error of the mean (SEM) is shown in figures, and P-values were calculated using a one-way or two-Comment citer ce document : Tronnet, S., Garcie, C., Brachmann, A. O., Piel, Oswald, E., Martin, P. (2017). High iron supply inhibits the synthesis of the genotoxin colibactin by pathogenic Escherichia coli through a non-canonical Fur/RyhB-mediated pathway. Pathogens and Disease, 75 (5). , DOI : 10.1093/femspd/ftx066 way ANOVA followed by a Bonferroni post-test unless otherwise stated. A P-value of less than 0.05 was considered statistically significant and is denoted by *. P<0.01 is denoted by ** and P<0.001 by ***, ns: not significant.

Results

High iron inhibits the production of colibactin genotoxin by the wild type strain

Recently we showed that the transcription of the clbA gene was modulated by iron (Tronnet et al. 2016). Therefore, we assessed the impact of iron supplementation on the production of the colibactin. We previously tested a panel of concentrations (1, 10, 100 and 1000 µM) (Tronnet et al. 2016). A dose effect of iron was observed on clbA expression. We chose to work with 100 µM of FeCl 3 , because with this concentration the maximum response of clbA transcription was reached without inducing toxicity in the bacteria (Tronnet et al. 2016).

The production of genotoxin colibactin was measured in wild type (WT) commensal strain M1/5 through bacteria-host cell interaction and subsequent quantification of megalocytosis (Figure 1A) and histone H2AX phosphorylation (Figure 1B) which both correlate with DNA double strand breaks resulting from the genotoxic effect of colibactin [START_REF] Nougayrède | Escherichia coli induces DNA double-strand breaks in eukaryotic cells[END_REF]. This revealed that the supplementation of the interaction medium (DMEM-Hepes) with FeCl 3 abrogated the genotoxic effect induced by the colibactin (Figure 1A,1B).

The mature colibactin genotoxin has not been fully characterized yet. However, a biosynthetic by-product derived from the colibactin assembly line, i.e. the N-myristoyl-D-asparagine moiety, has been recently characterized [START_REF] Bian | In vivo evidence for a prodrug activation mechanism during colibactin maturation[END_REF][START_REF] Brotherton | A prodrug resistance mechanism is involved in colibactin biosynthesis and cytotoxicity[END_REF]Vizcaino and Crawford 2015). This moiety is a prodrug motif that is cleaved during the late activation step (Dubois et al. 2011). We developed a LC-MS assay to quantitatively measure the amount of N-myristoyl-D-asparagine as a means to indirectly quantify the production of Comment citer ce document : Tronnet, S., Garcie, C., Brachmann, A. O., Piel, Oswald, E., Martin, P. (2017). High iron supply inhibits the synthesis of the genotoxin colibactin by pathogenic Escherichia coli through a non-canonical Fur/RyhB-mediated pathway. Pathogens and Disease, 75 (5). , DOI : 10.1093/femspd/ftx066 8 the genotoxin (Garcie et al. 2016) by E. coli M1/5 (Figure 1C). The results showed that in high iron conditions, the quantity of colibactin prodrug motif was abolished in the bacterial culture supernatant compared to the DMEM-Hepes (Figure 1C).

Altogether, these data showed that high iron inhibited the synthesis of colibactin in a wild type strain.

High iron inhibits both the transcription of clbA and the production of colibactin in Δfur mutant strains

Recently we showed that the transcription of the clbA gene and the production of colibactin was regulated by Fur. Therefore, we assessed the impact of high iron on the Fur-dependent regulation of clbA and colibactin.

E. coli strain Nissle 1917 harboring a transcriptional luciferase fusion with the clbA gene and its Δfur isogenic mutant were analyzed (Table 1, Homburg et al. 2007;Tronnet et al. 2016).

Both strain Nissle clbA-lux and strain Nissle clbA-lux Δfur were grown in chemically defined DMEM-Hepes and DMEM-Hepes supplemented with FeCl 3 (100μM). The transcription rate of clbA was determined as OD 600nm standardized-relative luminescence units (RLU) (Figure 2A). An increase of bioluminescence emission was observed during the bacterial growth, to reach a maximal value at early stationary phase of growth, followed by a decrease of transcription, for both strains (Figure 2A). In iron-supplemented conditions, the clbA transcription was strongly repressed in both strain Nissle clbA-lux and strain Nissle clbA-lux Δfur (Figure 2A).

The production of colibactin was measured in E. coli WT strains and Δfur mutant derivatives through bacteria-host cell interaction and subsequent quantification of megalocytosis in the absence and presence of iron supplementation (FeCl 3 , 100μM, Figure 2B). Two distinct genetic contexts of E. coli were investigated: Nissle 1917 (probiotic strain) and M1/5

Comment citer ce document : Tronnet, S., Garcie, C., Brachmann, A. O., Piel, Oswald, E., Martin, P. (2017). High iron supply inhibits the synthesis of the genotoxin colibactin by pathogenic Escherichia coli through a non-canonical Fur/RyhB-mediated pathway. Pathogens and Disease, 75 (5). , DOI : 10.1093/femspd/ftx066 9 (commensal strain). This revealed that in high iron interaction medium, the genotoxic effect induced by the colibactin was significantly diminished both in the WTs and the Δfur inactivated mutants.

Altogether, these data suggested that in Δfur mutants the iron-dependent clbA repression and the resulting decreased colibactin production could occur either through RyhB in a Furindependent manner or through an additional regulator.

High iron inhibits both the transcription of clbA and the production of colibactin in ΔryhB mutant strains

Recently we showed that the transcription of the clbA gene and the production of colibactin was regulated by RyhB. Therefore, we assessed the impact of high iron on the RyhBdependent regulation of clbA and colibactin.

E. coli strain Nissle clbA-lux and its ΔryhB isogenic mutant were analyzed (Table 1, Homburg et al. 2007;Tronnet et al. 2016). Both strain Nissle clbA-lux and strain Nissle clbA-lux ΔryhB were grown in DMEM-Hepes and DMEM-Hepes supplemented with FeCl 3 (100μM; Figure 3A). This revealed that in iron-supplemented conditions, the clbA transcription was repressed in both strains Nissle clbA-lux and Nissle clbA-lux ΔryhB (Figure 3A).

The production of colibactin was measured in E. coli WT strains and their ΔryhB mutant derivatives through bacteria-host cell interaction and subsequent quantification of megalocytosis in the absence and presence of iron supplementation (FeCl 3 , 100μM, Figure 3B). This revealed that in high iron interaction medium, the genotoxic effect induced by the colibactin was significantly diminished both in the WTs and the ΔryhB inactivated mutants.

Altogether, these data suggested that in a ΔryhB mutant the iron-dependent clbA repression and the resulting decreased colibactin production could occur through Fur in a ryhBindependent manner, or through an additional regulator.

Comment citer ce document : Tronnet, S., Garcie, C., Brachmann, A. O., Piel, Oswald, E., Martin, P. (2017). High iron supply inhibits the synthesis of the genotoxin colibactin by pathogenic Escherichia coli through a non-canonical Fur/RyhB-mediated pathway. Pathogens and Disease, 75 (5). , DOI : 10.1093/femspd/ftx066

10

A yet unknown factor is involved in the high iron-dependent regulation of colibactin production

To address the question of the potential involvement of an additional regulatory factor in the regulation of colibactin, ΔryhB Δfur double mutants were constructed in strains Nissle 1917 and M1/5, and were analyzed for their capacity to produce colibactin by quantification of megalocytosis (Figure 4A). This revealed that the quantity of colibactin produced in the double mutants was highly decreased, compared to the WT strains. The genotoxic effect measured in the double mutants was lower than that observed for the Δfur inactivated mutants, in each genetic context (Figure 4A).

To test the putative involvement of an additional regulatory factor in the high iron-dependent control of colibactin production, megalocytosis assay was performed with the double mutants ΔryhB Δfur in the presence of high iron (100μM, Figure 4B). This revealed that the colibactin-associated genotoxic effect was diminished by iron supplementation in the ΔryhB Δfur mutants.

Altogether, these findings suggested that an additional yet unknown regulatory factor was involved in the regulation of the colibactin in response to high iron conditions.

Discussion

Our previous study (Tronnet et al. 2016) highlighted that the mutation of ryhB led to an increase of colibactin production, whereas the mutation of fur resulted in a decrease of colibactin production in three different genetic contexts of E. coli belonging to the B2 phylogenetic group.

Here we showed that high iron decreases clbA transcription and colibactin production in wild type, ΔryhB, Δfur and ΔryhB Δfur strains. This suggests an iron-dependent synthesis of the genotoxin through a non-classic Fur/RyhB-mediated pathway. Nowadays, numerous reports Comment citer ce document : Tronnet, S., Garcie, C., Brachmann, A. O., Piel, Oswald, E., Martin, P. (2017). High iron supply inhibits the synthesis of the genotoxin colibactin by pathogenic Escherichia coli through a non-canonical Fur/RyhB-mediated pathway. Pathogens and Disease, 75 (5). , DOI : 10.1093/femspd/ftx066 11 support the four modes of Fur regulation, apo-and holo-Fur activation and repression, establishing a significant deviation from the classical model of Fur regulation [START_REF] Butcher | Structure and regulon of Campylobacter jejuni ferric uptake regulator Fur define apo-Fur regulation[END_REF]). In the same manner, ryhB was reported to be regulated in a Fur-dependent and Furindependent manner, and to repress or activate genes [START_REF] Massé | Effect of RyhB small RNA on global iron use in Escherichia coli[END_REF][START_REF] Papenfort | Target activation by regulatory RNAs in bacteria[END_REF].

Colibactin is a virulence factor involved in neonatal systemic infections [START_REF] Mccarthy | The genotoxin Colibactin is a determinant of virulence in Escherichia coli K1 experimental neonatal systemic infection[END_REF], and is also associated with the development of colorectal cancer [START_REF] Arthur | Intestinal inflammation targets cancerinducing activity of the microbiota[END_REF][START_REF] Cougnoux | Bacterial genotoxin colibactin promotes colon tumour growth by inducing a senescence-associated secretory phenotype[END_REF][START_REF] Cougnoux | Bacterial genotoxin colibactin promotes colon tumour growth by inducing a senescence-associated secretory phenotype[END_REF]. Our work could make a link between iron concentration in the gut/the blood and E. coli-mediated carcinogenesis/systemic infections. Given the existence of a gradient of iron concentration from the lumen to the intestinal epithelial cell, it is conceivable that the fine-tuning of clbA expression allows the production of colibactin when the pathogenic E. coli is located in an appropriate site in the gut. This iron-dependent regulation could also be a reason for the emergence of B2 group E. coli, as we know there is a crosstalk between colibactin and siderophores [START_REF] Martin | Interplay between siderophores and colibactin genotoxin biosynthetic pathways in Escherichia coli[END_REF].

Integrating the regulation of virulence factors, such as siderophores and colibactin, into networks that respond to specific environmental signals, such as the local iron concentration and the balance between Fur and RyhB could result in an accurate production of colibactin and siderophores, so that the bacteria can adapt to the competitive environment that is the gut, or in the blood.

This work highlights the complex regulatory mechanism of an important bacterial virulence and carcinogenesis factor by major regulators of bacterial iron homeostasis, and provides insights into the consequences associated with high iron supply (graphical abstract). Further understanding of how E. coli senses environments to regulate the biosynthesis of colibactin may uncover novel pathways for the development of potential targets against this pathogenic enterobacteria. A. The production of colibactin by E. coli strain M1/5 was determined by quantification of megalocytosis as previously described [START_REF] Martin | Interplay between siderophores and colibactin genotoxin biosynthetic pathways in Escherichia coli[END_REF]. HeLa cells were infected with E. coli, in the presence of 100µM FeCl 3 or not, then washed 4h post infection. After infection, the cells were incubated for 72h with gentamicin before protein staining with methylene blue.

Figure 1 .

 1 Figure 1. High iron inhibits the production of colibactin genotoxin by the wild type

  The staining was quantified by acid-extraction of the methylene blue and measurement of optical density at 660nm (OD 660 nm ). Multiplicity of infection (number of bacteria per cells): MOI = 200 and 100. B. The production of colibactin by E. coli strain M1/5 was determined by quantification of H2AX phosphorylation. HeLa cells were infected with E. coli, in the presence of 100µM FeCl 3 or not, then washed 4h post infection. The cells were then incubated 3h in DMEM medium supplemented with gentamicin before fixation, then Comment citer ce document : Tronnet, S., Garcie, C., Brachmann, A. O., Piel, Oswald, E., Martin, P. (2017). High iron supply inhibits the synthesis of the genotoxin colibactin by pathogenic Escherichia coli through a non-canonical Fur/RyhB-mediated pathway. Pathogens and Disease, 75 (5). , DOI : 10.1093/femspd/ftx066 18 permeabilized and labeled for DNA and phosphorylated H2AX histone (γ-H2AX) using an In-Cell Western method (Martin et al. 2013). MOI = 50. C. The colibactin prodrug motif Nmyristoyl-D-asparagine produced by E. coli strain M1/5 was quantified by LC-MS. Bacteria were cultivated at 37°C for 18h in DMEM-Hepes or DMEM-Hepes supplemented with FeCl 3 (100µM), and N-myristoyl-D-asparagine was quantified in culture supernatants and pellets by LC-MS using N-myristoyl-l-asparagine (isomer of the N-myristoyl-D-asparagine colibactin prodrug motif) as a standard. The results were normalized to the bacterial biomass and are presented as quantity of N-myristoyl-D-asparagine (ng/ml). Data represented in the graph were obtained from three biological replicates and two independent experiments. NI: not infected. Statistical analysis: unpaired t-test. ***: p < 0.001, **: p < 0.01. ns: not significant.

Figure 2 .

 2 Figure 2. High iron inhibits the transcription of clbA and the production of colibactin

Figure 3 .

 3 Figure 3. High iron inhibits the transcription of clbA and the production of colibactin

Figure 4 .

 4 Figure 4. A yet unknown factor is involved in the high iron-dependent regulation of

Table 1 .

 1 Strains and primers used in this study. Nissle clbA-lux Luciferase fusion of the gene clbA, Kan R Homburg et al. 2007 Nissle clbA-lux ΔryhB ryhB mutant of strain Nissle clbA-lux, Kan R Cm R Tronnet et al. 2016 Nissle clbA-lux Δfur fur mutant of strain Nissle clbA-lux, Kan R Cm R Tronnet et al. 2016 Nissle clbA-lux ΔryhB Δfur ryhB fur double mutant of strain Nissle clbA-lux Kan R Cm R

	E. coli strain

Funding sources

This work was supported by the Agence Nationale de la Recherche (France) [grant ANR-13-BSV3-0015-02, ANR-13-BSV1-0028-01].

Conflict of interest

The authors have no conflict of interest to declare.