Pierre Hohmann
email: pierre.hohmann@fibl.org

Benedikt Haug

Monika M Messmer

Jérôme Enjalbert

Isabelle Goldringer

Emma Forst

Timothée Flutre

Tristan Mary-Huard

Advances in Breeding for Mixed Cropping -Incomplete Factorials and the Producer/Associate Concept

Keywords: mixed cropping, intercropping, breeding, general mixing ability, producer/associate concept, incomplete factorial design, biological interaction, simulations

Mixed cropping has been suggested as a resource-efficient approach to meet high produce demands while maintaining biodiversity and minimizing environmental impact. Current breeding programs do not select for enhanced general mixing ability (GMA) and neglect biological interactions within species mixtures. Clear concepts and efficient experimental designs, adapted to breeding for mixed cropping and encoded into appropriate statistical models, are lacking. Thus, a model framework for GMA and SMA (specific mixing ability) was established. Results of a simulation study showed that an incomplete factorial design combines advantages of two commonly used full factorials, and enables to estimate GMA, SMA, and their variances in a resource-efficient way. This model was extended to the Producer (Pr) and Associate (As) concept to exploit additional information based on fraction yields. It was shown that the Pr/As concept allows to characterize genotypes for their contribution to total mixture yield, and, when relating to plant traits, allows to describe biological interaction functions (BIF) in a mixed crop. Incomplete factorial designs show the potential to drastically improve genetic gain by testing an increased number of genotypes using the same amount of resources. The Pr/As concept can further be employed to maximize GMA in an informed and efficient way. The BIF of a trait can be used to optimize species ratios at harvest as well as to extend our understanding of competitive and facilitative interactions in a mixed plant community. This study provides an integrative methodological framework to promote breeding for mixed cropping.

INTRODUCTION

Climate change, such as rising global temperatures and climatic volatility are predicted to jeopardize future agricultural productivity [START_REF] Rahmstorf | Increase of extreme events in a warming world[END_REF]. The current strategies to produce stable and high yields, e.g., by the application of mineral fertilizer, are of limited future use since they themselves are a contributor to these changing climatic parameters [START_REF] Thompson | Acceleration of global N 2 O emissions seen from two decades of atmospheric inversion[END_REF]. Thus, alternative approaches to achieve high and stable yields while maintaining biodiversity and minimizing environmental impact have to be developed. Mixed cropping is the simultaneous cultivation of two or more crops on the same field. Especially legume/non-legume species mixtures have been proposed to achieve a higher per area production and profitability and higher yield stability with less or no external inputs [START_REF] Bedoussac | Ecological principles underlying the increase of productivity achieved by cereal-grain legume intercrops in organic farming[END_REF][START_REF] Raseduzzaman | Does intercropping enhance yield stability in arable crop production? A meta-analysis[END_REF][START_REF] Wendling | Specific interactions leading to transgressive overyielding in cover crop mixtures[END_REF][START_REF] Viguier | Yield gap analysis extended to marketable grain reveals the profitability of organic lentil-spring wheat intercrops[END_REF]. Good pairs of complementary species have already been identified, such as combinations of corn (Zea mays L.) with cowpea (Vigna unguiculata L., [START_REF] Ofori | Maize/cowpea intercrop system: effect of nitrogen fertilizer on productivity and efficiency[END_REF], with common bean (Phaseolus vulgaris L., [START_REF] Hoppe | Entwicklung von Energiemaissorten für die Mischkultur mit Stangenbohnen[END_REF][START_REF] Starke | Selektion von Stangenbohnensorten (Phaseolus vulgaris L.) für den Mischanbau mit Mais[END_REF], and with faba bean (Vicia faba L., [START_REF] Li | Syndromes of production in intercropping impact yield gains[END_REF], as well as small grain cereals such as barley (Hordeum vulgare L., Hauggaard-Nielsen et al., 2001) with pea (Pisum sativum L.) or wheat (Triticum aestivum L.) with faba bean [START_REF] Agegnehu | Yield performance and land-use efficiency of barley and faba bean mixed cropping in Ethiopian highlands[END_REF] or with lentil (Lens culinaris MEDIK., [START_REF] Viguier | Yield gap analysis extended to marketable grain reveals the profitability of organic lentil-spring wheat intercrops[END_REF].

The choice of the genotypes best suited to mixing within each species is not straight-forward, and a robust strategy for evaluating and identifying the right mixing partners from among a large number of candidates is essential to selecting the components of mixtures and improving the mixing ability of each species. Selection efficiency for mixed cropping yield under pure stand has been reported to be moderate or low highlighting the value dedicated breeding efforts for mixed cropping [START_REF] De Oliveira Zimmermann | Breeding for yield, in mixtures of common beans (Phaseolus vulgaris L.) and maize (Zea mays L.)[END_REF][START_REF] O'leary | Uncovering corn adaptation to intercrop with bean by selecting for system yield in the intercrop environment[END_REF][START_REF] Annicchiarico | Do we need specific breeding for legumebased mixtures?[END_REF]. Well performing genotypes should display a high general mixing ability (GMA), i.e., lead to a high total mixture yield performance across several potential mixing partners, and a low variance in specific mixing ability (SMA), i.e., little or no specific interaction with individual mixing partners in that respect. In order to develop efficient breeding strategies for crop mixtures of two species, trial designs must be developed that allow precise estimation of GMA and SMA variance.

Current trial designs apply a factorial setup, combining m genotypes of species one and n genotypes of species two. With increasing numbers of genotypes of both mixing partners, factorial designs quickly result in an unfeasibly high number of experimental plots. Therefore, often specific dimensions of both crop species are used: depending on the question to be addressed either (i) factorials of equal (or similar) dimensions of m and n with a small to medium number of genotypes [START_REF] Annicchiarico | Interference effects in white clover genotypes grown as pure stands and binary mixtures with different grass species and varieties[END_REF]Hauggaard-Nielsen and Jensen, 2001) or (ii) factorials of different dimensions for m and n [START_REF] Annicchiarico | Breeding white clover for increased ability to compete with associated grasses[END_REF][START_REF] Hoppe | Entwicklung von Energiemaissorten für die Mischkultur mit Stangenbohnen[END_REF][START_REF] Starke | Selektion von Stangenbohnensorten (Phaseolus vulgaris L.) für den Mischanbau mit Mais[END_REF] are employed. The former allows GMA and SMA estimations of both crop species involved, the latter emphasizes one species over the other and is comparable to the topcross designs, used in hybrid breeding. With advances in mixed modeling statistical software (such as with the GNU R packages "lme4" or "SOMMER"), analyses of largely incomplete datasets are possible [START_REF] Bates | Fitting linear mixed-effects models using lme4[END_REF][START_REF] Covarrubias-Pazaran | Genome-Assisted Prediction of Quantitative Traits Using the R Package sommer[END_REF][START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]. Incomplete factorial designs have been suggested to mitigate the limitations of (i) and (ii) by expanding the numbers of m and n while maintaining a feasible number of experimental plots. Previously, they have been applied to assess GMA and SMA effects in wheat (Triticum aestivum L.) cultivar mixtures [START_REF] Forst | A generalized statistical framework to assess mixing ability from incomplete mixing designs using binary or higher order variety mixtures and application to wheat[END_REF] and found recent application in genomic prediction in corn (Zea mays L.) hybrid breeding [START_REF] Seye | Revisiting hybrid breeding designs using genomic predictions: simulations highlight the superiority of incomplete factorials between segregating families over topcross designs[END_REF].

Another important area to advance breeding for mixed cropping focuses on exploiting information that is contained in the fraction yields of mixed crops via the application of the producer (Pr) and associate (As) concept [START_REF] Wright | Selection for improved yield in inter-specific mixtures or intercrops[END_REF][START_REF] Goldringer | Adjustment for competition between genotypes in single-row-plot trials of winter wheat (Triticum aestivum)[END_REF][START_REF] Annicchiarico | Do we need specific breeding for legumebased mixtures?[END_REF]. In the mixed cropping context, the Pr effect, sometimes also referred to as direct effect, is the capacity of a genotype to influence its own yield in a mixture, while the As effect is its capacity to influence the yield of its companion crop or variety [START_REF] Annicchiarico | Do we need specific breeding for legumebased mixtures?[END_REF]. As laid out by [START_REF] Wright | Selection for improved yield in inter-specific mixtures or intercrops[END_REF], [START_REF] Forst | Développement de Méthodes D'estimation de L'aptitude au Mélange Pour la Prédiction des Performances et la Sélection de Mélanges Variétaux Chez le Blé Tendre et co-Conception D'idéotypes De Mélanges Adaptés À L'agriculture Biologique[END_REF][START_REF] Sampoux | Which recurrent selection scheme to improve mixtures of crop species? theoretical expectations[END_REF], the Pr and As effects of a given genotype sum up to its GMA effect. It has been applied to single row experiments in breeding nurseries [START_REF] Goldringer | Adjustment for competition between genotypes in single-row-plot trials of winter wheat (Triticum aestivum)[END_REF] or wheat cultivar mixtures [START_REF] Forst | Développement de Méthodes D'estimation de L'aptitude au Mélange Pour la Prédiction des Performances et la Sélection de Mélanges Variétaux Chez le Blé Tendre et co-Conception D'idéotypes De Mélanges Adaptés À L'agriculture Biologique[END_REF], and to the mixed cropping context [START_REF] Wright | Selection for improved yield in inter-specific mixtures or intercrops[END_REF][START_REF] Sampoux | Which recurrent selection scheme to improve mixtures of crop species? theoretical expectations[END_REF]. Separated yield data enables either uni-or bivariate (or multivariate) analysis, i.e., joint analysis of the two (or multiple) fraction yields. In clinical psychology as well as in livestock breeding, multivariate analysis procedures have successfully been applied in situations where traits were correlated, e.g., due to pleiotropy, and yielded higher precision for QTL detection than univariate approaches [START_REF] Sørensen | A comparison of bivariate and univariate QTL mapping in livestock populations[END_REF][START_REF] Meier | Joint analysis of psychiatric disorders increases accuracy of risk prediction for schizophrenia, bipolar disorder, and major depressive disorder[END_REF]. In mixed cropping conditions, it can be assumed that errors of measurements are generally negatively correlated between the two crops, e.g., via compensation effects. Often in mixed cropping, one of the species is at a competitive disadvantage and genotypes of the species that is very non-competitive generally express a low GMA in mixtures [START_REF] Annicchiarico | Interference effects in white clover genotypes grown as pure stands and binary mixtures with different grass species and varieties[END_REF][START_REF] Corre-Hellou | Interspecific competition for soil N and its interaction with N2 fixation, leaf expansion and crop growth in pea-barley intercrops[END_REF]. However, the competitive ability of genotypes is obscured when only whole mixture yield is observed [START_REF] Annicchiarico | Do we need specific breeding for legumebased mixtures?[END_REF]. The assessment of fraction yields is not only important to identify these competitive abilities, but can also be applied to optimize a mixture toward a specific ratio, e.g., for feed nutrition or legume subsidies reasons, as is for example the case in Switzerland (Bundesamt für Landwirtschaft BLW, 2019). Thus, shaping community performance and composition via traits can be of interest in the breeding process. Furthermore, high genetic correlations between certain traits and mixture yield would allow indirect selection based on most important key traits for Pr and As effects.

The choice of an efficient trial design, the choice of an efficient analysis method and the assessment of yield proportions are interrelated topics. They provide the potential to be combined in an integrated approach to promote breeding for mixed cropping. While some published work focuses on the parallel genetic improvement of two species [START_REF] Sampoux | Which recurrent selection scheme to improve mixtures of crop species? theoretical expectations[END_REF], many publications rely on the improvement of one species at a time and do not take the potential of analyzing separated yield data into account. Thus, the objectives of this study were to (i) develop a model to estimate GMA and SMA variances of binary species mixtures, and to compare different experimental designs for their usefulness in estimating these parameters, (ii) subdivide GMA into Pr and As effects in order to categorize cultivars' influence on mixture yield and to compare the precision of a uni-versus a bivariate approach in estimating Pr, As, and error (co)variances, and (iii) establish a concept to link plant traits to biological interactions between involved species in mixture.

MATERIALS AND METHODS

To exemplify the case of a mixed crop, a hypothetic binary mixture of a legume species (pea) and a non-legume species (barley) will be used in this study.

The GMA Model of Total Mixture Yield

Mixture yield can be expressed with the following model:

y ijk = µ + r k + G p i + G b j + S ij + E ijk (1)
with y ijk the total mixture yield of the i-th pea cultivar mixed with the j-th barley cultivar in the k-th block, µ the intercept of mixture yields, r k the effect of the k-th block (replication), G p i and G b j the GMA effects the i-th pea cultivar and the j-th barley cultivar, respectively, S ij the SMA effect, i.e., interaction, of the i-th pea cultivar with the j-th barley cultivar and E ijk the error term. This model-framework was used to compare four different trial designs (A-D), comprising three full ("f ") factorials, with all possible pairwise combinations present, and an incomplete ("i") factorial with only a subset of all possible pairwise combinations present (Figure 1). In the following, m is the number of barley cultivars and n the number of pea cultivars, used in the design. Design A, using 240 plots per replicate (m = 8 and n = 30), is the most resource-expensive design. Designs B, C, and D were using only approximately 25% of the resources of A with, respectively, 64 (m = 8, n = 8), 60 (m = 8, n = 30), and 60 (m = 2, n = 30) plots per replicate, while sharing commonalities with design A. Designs B and C are both also full factorials (with equal and unequal dimensions of m and n) and design D shares the same size of m and n with design A while being an incomplete factorial. Design C (2 × 30f) full factorial has similarities with a top cross design, used in early stages of a hybrid breeding program. Design D (8 × 30i) was constructed using four independently randomized 8 × 8 latin squares of which only two "entries" were used for the mixtures. This ensured that (a) every barley was combined with eight pea cultivars, and every pea was combined with two different barley cultivars and (b) confounding, i.e., two or more pea cultivars sharing the same two barley cultivars was minimized.

For the comparison of the four designs, datasets with total mixture yield-data ("total yield setting") were simulated with the following models for a "SMA present" and a "SMA absent" simulation, respectively. For the SMA absent simulation, S ij was set to zero in the model in formula (1). Simulations were performed according to the following procedure, using parameter settings in the same order of magnitude as empirical values from preliminary trials (Haug et al., in preparation) to produce data that is as close to empirical data as possible. The intercept of mixture yield was set to 38.8 dt/ha. For each of the following parameters, the corresponding effects were drawn from their respective probability distributions: block effects r k were drawn from a normal distribution with a mean of 0 and a variance of 2, i.e., N (0, 2). Pea GMA effects G p i were drawn from N (3, 0), barley GMA effects G b j were drawn from N (0, 5) and SMA effects S ij were drawn from N (0, 5) for the SMA present simulation and from N (0, 0) for the SMA absent simulation, i.e., the effect size was set to zero. Errors e ijk were drawn from N (0, 5). For each simulation run, effects were drawn anew. "True" values (i.e., values used to simulate data) of each effect were saved after each simulation run for later comparison with the estimated parameter values, e.g., the true GMA effects of the pea cultivars, were stored for later comparison with GMA effects estimated by the best linear unbiased predictors (BLUPs) received by the mixed model to analyze the simulated data. Design D normally would have 64 pairwise combinations. The realistic case of missing/unusable genotypes was assumed, by excluding two pea lines in our case, thus resulting in 60 pea-barley combinations. For each of the four trial designs n = 1,000 data sets were simulated, for the SMA present, as well as the SMA absent simulation, i.e., 8 × 1,000 simulated data-sets in total. Each dataset comprised two replicates (blocks). All simulations and subsequent analyses were done using GNU R (R Core Team, 2019). The R-code used for simulation and analysis is publically available (Haug, 2020a).

The Pr/As Model of Fraction Yield of Each Species

Effects on separated yield data, i.e., pea and barley fraction yields, can be described by a model, containing Pr and As effects:

y p ijk = µ p + r p k + P p i + A b j + E p ijk
(2)

y b ijk = µ b + r b k + P b j + A p i + E b ijk (3)
with y p ijk the fraction yield of pea from the combination of the i-th pea cultivar with the j-th barley cultivar in the k-th block, µ p the intercept of pea fraction yields, P p i the effects of the k-th block on pea fraction yield, P p i the Pr effects and A b j the As effects of the i-th pea and the j-th barley cultivar, respectively, and E p ijk the error for the fraction yield of the i-th pea with the j-th barley in the k-th block. Each dataset comprised two replicates (blocks). Parameters apply in analogy for barley fraction yields in formula 3. Interactions between Pr and As effects are ignored. Since total mixture yield y ijk decomposes into (y p ijk + y b ijk ), also the other parameters can be decomposed:

µ into µ p + µ b , r k into (r p k + r b k ), G p into (P p + A p ), G b into (P b + A b )
, and E ijk into (E p ijk + E b ijk ). Hence, formula 1 (without S ij can be rewritten as in formula 4.

y p ijk + y b ijk = µ p + µ b + r p k + r b k + P p i + A p i + P b j + A b j +E p ijk + E b ijk (4)
P p i A p i ∼ N 2 0 0 , 6.1 -4.1 -4.1 7.5
(5)

P b j A b j ∼ N 2 0 0 , 6.1 -2.7 -2.7 1.5 (6) P b j A b j ∼ N 2 0 0 , 6.1 3.6 3.6 8.4 (7)
Decomposition of GMA into Pr and As effects is illustrated in Figure 2. Separated yield data ("fraction yield setting") were simulated to compare the precision of a univariate versus a bivariate analysis approach for design D (Figure 1). As in the total yield setting, two blocks per data set were simulated. In contrast to the total yield setting, for each "plot" the separated yields of pea and barley were simulated with a mean pea yield of 18.3 dt/ha and a mean barely yield of 19.1 dt/ha, according to previously mentioned preliminary experimental data. The following settings were used: block effects r p k and r b k were both drawn from N (0, 2) as in the total yield setting. Pea Pr effects P p i were drawn from N (0, 6.1) and As effects A b j from N (0, 7.5). Barley Pr effects P b j were drawn from N (0, 6.1) and As effects A p i from N (0, 1.5), errors E p ijk and E b ijk were drawn for pea yield from N (0, 8.4) and for barley yield from N (0, 6.1). Correlations of pea Pr effects with pea As effects were set to -0.61, and correlation of barley Pr with barley As effects were set to -0.81. Within-plot error-correlations, i.e., the correlation of errors of barley yields with the errors of the corresponding pea yields in the same plot, were set to 0, -0.5, and -0.9, respectively, to create three different error-correlation scenarios. Withinplot error-correlations and Pr/As correlations were translated into co-variances. The variance-covariance matrices of Pr/As effects of pea, barley and the errors are shown in formulas 5, 6, and 7. Pr and As effects as well as errors were drawn from a distribution that follows the law of a multivariate normal distribution, using the function "mvrnorm" from the R-package "MASS" [START_REF] Venables | Modern Applied Statistics with S[END_REF], and the covariances shown in formulas 5-7. In total, three data sets for each of the three different error correlation settings with 1,000 simulations were created. Simulated data from the fraction yield setting was analyzed using (i) a univariate approach with models equal to those used to simulate the data (formulas 2, 3) and (ii) a bivariate approach (formula 4) in which the two dependent variables were analyzed jointly (Covarrubias-Pazaran, 2018). In addition to the parameters estimated by the univariate approach, the bivariate approach also estimates the before mentioned covariances. Both approaches were done as mixed models where block-effects were considered as fixed and all other effects as random, assuming independent and identically distributed random variables. The uni-and the bivariate analyses were done with the "mmer" function of the R-package "SOMMER" [START_REF] Covarrubias-Pazaran | Genome-Assisted Prediction of Quantitative Traits Using the R Package sommer[END_REF], 2018). Estimates of the model parameters, e.g., estimated GMA or Pr variances of pea, and BLUPs for the genetic effects, e.g., BLUPs of GMA effects, were saved for later analysis for each of the 1,000 datasets per setting. Depending on the analysis approach, each dataset yielded a different set of BLUPs for the pea and barley Pr and associate effects. For n = 1,000 analyses, Pearson correlations between the two sets of BLUPs and the true value were computed, Fisher-ztransformed, averaged for each coefficient and transformed back. T-tests between the mean correlation of the univariate and the bivariate approach within each parameter were conducted to compare the approaches for their accuracy to estimate the true effect values. The R-code used for simulation and analysis is publically available [START_REF] Haug | Simulation code for comparison of uni-and bivariate analysis in mixed cropping experiments[END_REF].

Trait versus GMA/Pr/As Analyses for the Characterization of Biological Interaction Functions (BIFs)

Beyond the purely statistical treatment of the data described above, the relationships between a fictive explanatory trait and the GMA/Pr/As variables were investigated. This explanatory trait was set into relation with the GMA, Pr, and As effects of equally fictive genotypes. Nine possible scenarios of trait-GMA, trait-Pr effect, and trait As-effect relationships were investigated, allowing the categorization of traits according to their biological interaction function in a mixed-crop plant community. Exemplary scatter plots with n = 100 simulated genotypes and positive, null and negative relationships between trait and GMA, trait and Pr, and trait and As effects were created to suggest a simple visual analysis of these relationships. Potential symbiotic trait functions were associated to the corresponding functions relationships.

RESULTS

Incomplete Design Yields Comparable Estimates to Full Factorial Designs

Four experimental designs were compared for their ability to estimate GMA and SMA variances as well as estimating genotypic effects (BLUPs) correctly in two different simulations, a "SMA-present" and a "SMA-absent" simulation.

Over both scenarios, the precision of the estimates increased with experimental resource input (Table 1). With design A (8 × 30f), utilizing 240 experimental plots per replicate, the narrowest CIs among the four designs were received, more narrow than the ones of designs B (8 × 8f), C (2 × 8f), and D (8 × 30i), which were using only 64, 60, and 60 experimental plots per replicate, respectively (Table 1). Among the latter, only minor differences in CIs were observed (except lower reliability on GMA variance of barley in design C and on GMA variance of pea in design B in the SMA-absent simulation). Besides GMA variance of barley of design C in the SMA-present simulation and the GMA variance of pea of design B in the SMA-absent simulation, certain parameters were estimated similarly well with designs B, C, and D as with the benchmark design A. Barley GMA variance of designs B and D showed similar CIs compared to design A for both the SMA-present and the SMA-absent simulation, whereas design C, which only uses two instead of eight barley cultivars, estimated barley GMA variance less precisely (CI of ±0.43 in both scenarios). In addition, pea GMA variance of designs C and D of the SMA-absent simulation were similarly precisely estimated compared with design A, whereas design B estimated this parameter with lower precision (CI of ±0.11).

Besides the variation of estimates, a check of the correct estimation of the size of the parameter itself revealed that for the SMA-present simulation (Table 1), all four experimental designs accurately estimated GMA, SMA, and error variances, with all confidence intervals (CIs) of means overlapping the true values, except for the SMA variance of design C.

For the SMA-absent simulation, pea and barley GMA variances were mostly accurately estimated in the four designs with significant but small underestimations of pea GMA variance by designs C and D. SMA variances in this simulation were significantly overestimated and error variances were significantly underestimated for all four designs. Compared with the mean SMA variance of the benchmark design A of 0.14, in this The four designs employ different amounts of resources in terms of experimental plots per replicate, with design A being the benchmark. Numbers after ± denote the 95% confidence intervals for the mean of 1,000 variance estimates. Different letters following correlation values between best linear unbiased predictors (BLUPs) and true effects indicate significant differences (Tukey-HSD; p < 0.05).

1 the different designs as described in simulation designs B, C, and D showed significantly higher SMA variances with 0.27, 0.35, and 0.36, respectively. Similarly, the mean error variance of design A was with 4.88 significantly higher (and thus closer to the truth of 5.0) than the error variances of designs B, C, and D, with means between 4.76 and 4.80. When comparing the four designs for their correlation of BLUPs with the truth value, (Table 1) the correlation coefficients of benchmark design A were among the highest across all estimated BLUPs. However, correlations of BLUPs of pea GMA effects with the true effects of design B were similar compared with the benchmark design A for both simulations. For the correlations of BLUPs for pea GMA with their true values of designs C and D showed significantly lower correlation coefficients compared with design A and B for both SMA simulation models. However, this difference was less apparent in the SMA-absent simulation with correlation coefficients of 0.84 and 0.83 of designs C and D compared with 0.95 and 0.96 of designs A and B, respectively. All correlation coefficients differed significantly from each other (p < 0.001). For barley GMA, design B and D showed with 0.93 and 0.91 high mean correlation coefficients that were similar to the mean of design (0.98). Correlation coefficients for SMA effects were in the range of the correlation coefficients for pea GMA effects and barley GMA effects, with values between 0.67 (design D) and 0.77 (design A).

The Pr/As Concept Allows to Characterize Cultivars' Contribution to Mixture Yield

In Figure 3, thirty simulated cultivars with their Pr and As effects are shown. Pr effects range from -5.0 to +7.0 and As effects range from -5.4 to +4.3. As effect has to be read as the effect of the species (e.g., pea) on the yield (or any other trait) of its companion species (e.g., barley). Since the Pr and As effects of a cultivar sum up to its GMA, cultivars that lie on the line with slope -1 and intercept 0 have a GMA of zero, those above this line have a positive and below a negative GMA. The cultivars thus can be grouped into six sectors, U, V, W, X, Y, and Z, with U-W having positive GMA due to a high Pr effect that offsets a negative As effect (sector U), both a positive Pr and As effect (sector V) and a positive As effect that makes up for a negative Pr effect (sector W). On the other hand, cultivars below the identity line have a negative GMA with a positive Pr effect that does not compensate for a negative As effect (sector X), both negative Pr and As effects (sector Y) and a negative Pr effect which is not offset by their positive As effect (sector Z). These six sectors allow to characterize and differentiate the mixing ability of the pea cultivars.

Comparing the uni-and bivariate approach using the data from the fraction yield setting, correlation of errors of 0, -0.5, and -0.9 resulted in very similar parameter estimates (Supplementary Table 1). The "error correlation of -0.5 scenario" was used to analyze separated fraction yield data of both models in more detail. Both analysis approaches, uniand bivariate, produced unbiased results of parameter estimates, i.e., all 95% confidence intervals of the estimates contained the true values for both approaches (see Supplementary Table 2).

Estimates did not differ significantly between the univariate and the bivariate analysis approach. CIs, used as a measure for precision, differed only by 0.01 or not at all. However, only the bivariate model allows to estimate the correlation between Pr and As effects.

Pr/As-Trait Relationships Uncover Biological Interaction Functions (BIFs) of Traits

The GMA, Pr, and As effects on total or fraction yield do not reveal the underlying biological processes or traits that influence the mixing ability. Yet, the examination of relationships between a fictive explanatory trait and Pr/As effects on fraction yield fills this lack by defining nine potential biological interaction functions (BIFs) of a given trait that underlie the GMA-trait pattern (Figure 3). This GMA-trait relationship can be positive (+), absent (0) or negative (-). However, the GMA-trait relationship is subdivided in its underpinned three possible Pr-trait/As-trait relationships. These can then be interpreted in terms of BIFs: commensalism (Pr+/As0, Pr0/As+, i.e., trait FIGURE 3 | Best linear unbiased predictors of producer (Pr) and associate (As) effects of 30 pea cultivars of a bivariate analysis of simulated data in which an 8 barleys × 30 peas incomplete factorial design was used (Design C). Pr and As effects represent the yield effects in dt/ha of cultivars of a focal species (here pea) on partial yields of itself (Pr) or on the associated species (As; here barley; read "As effect of pea on barley yields"). Data taken from a randomly chosen simulated data set of the Pr/As data. The sum of the Pr and As effects of a cultivar equals its GMA effect, thus, the line with slope -1 and intercept 0 separates genotypes with positive (above) and negative (below) GMA. Genotypes can have a positive GMA by either a high Pr effect that offsets for a negative. As effect (sector U), both positive Pr and As effects (sector V) or a high As effect that offsets for a negative Pr effect (sector W; consequently opposite for sectors X, Y, and Z).

will profit only one species), mutualism (Pr+/As+, i.e., trait promotes both species), antagonism (Pr+/As-, Pr-/As+, i.e., trait promotes one species but hampers second species), neutralism (Pr0/As0, i.e., trait does not affect any of the two species), amensalism (Pr0/As-, Pr-/As0, i.e., trait is hampering only one species), and competition (Pr-/As-, i.e., trait is hampering both species). This more detailed correlations will allow to identify key traits that are important for good mixing ability for a crop and can contribute to indirect selection.

DISCUSSION

The overall goal of this study was to develop a novel framework for breeding for mixed cropping by (i) formulating models for mixed cropping, suggesting experimental designs and analysis methods for mixed cropping experiments, (ii) proposing extensions to the use of the Pr/As concept, and (iii) linking the latter to traits in order to uncover the biological interaction function (BIF) of traits.

Incomplete Designs to Increase Selection Intensities

With all four designs, GMA, SMA, and error variances were overall correctly estimated, i.e., with little or no bias. As expected, the low-resource designs B, C, and D showed a slightly lower precision (i.e., higher CIs). The comparison of the three lowresource designs, incomplete factorial design (D) versus the two full factorial designs B and C, revealed similar estimations of GMA, SMA, and error variances. Design C (with two barley genotypes) does not allow a meaningful estimation of GMA variance. Thus, for an estimation of both species' GMA variance and SMA variance, designs B or D are preferred, which is in line with the suggestion of [START_REF] Annicchiarico | Do we need specific breeding for legumebased mixtures?[END_REF]. [START_REF] Seye | Revisiting hybrid breeding designs using genomic predictions: simulations highlight the superiority of incomplete factorials between segregating families over topcross designs[END_REF] used an incomplete factorial design for hybrid testing that was created by crossing one inbred line of one pool with one inbred line from the opposite other pool and compared to a classical topcross design where all inbred lines of pool 1 were crossed with the same line of pool 2. In the incomplete design, estimates for general combining ability (GCA) for both parents of a hybrid cannot be disentangled and are identical. Nonetheless, they emphasize, that if one only considers selection among the tested lines of the training set, the incomplete factorial (similar to design D of this study) always outperformed a topcross design (similar to design C) in terms of genetic gain, since twice the amount of genotypes (similar to barley in our case) could be tested and consequently selection intensity could be twice as high.

Correlations of true with estimated pea GMA values are lower in the incomplete design D compared with design B. In the absence of SMA, however, correlations come close to those of design B and even high-resource design A. This suggests that incomplete factorials can combine the advantages and applications of both designs with equal or similar dimensions of m and n and of designs with unequal dimensions of m and n, the latter being similar to topcross designs in hybrid breeding. At similar resource requirements, incomplete factorials allow more genotypes to be tested without a substantial loss of GMA precision and accuracy. This will allow to increase selection intensity similar to the example shown in maize hybrid breeding [START_REF] Seye | Revisiting hybrid breeding designs using genomic predictions: simulations highlight the superiority of incomplete factorials between segregating families over topcross designs[END_REF]. Moreover, testing larger sets of genotypes allow to exploit larger genetic variance of a given species. Since selection gain depends on both intensity and genetic variance, incomplete designs have a great potential to increase selection gain when breeding for mixed cropping.

Incomplete Designs for Early and Later Stages of Breeding for Mixed Cropping

Besides having been suggested for calibration of genomic prediction models for hybrid breeding [START_REF] Seye | Revisiting hybrid breeding designs using genomic predictions: simulations highlight the superiority of incomplete factorials between segregating families over topcross designs[END_REF], incomplete designs have been used to estimate GMA and SMA effects in wheat cultivar mixtures [START_REF] Forst | A generalized statistical framework to assess mixing ability from incomplete mixing designs using binary or higher order variety mixtures and application to wheat[END_REF]. The findings of [START_REF] Forst | A generalized statistical framework to assess mixing ability from incomplete mixing designs using binary or higher order variety mixtures and application to wheat[END_REF] could be applied to hybrid breeding as well as to breeding for mixed cropping: in the early development of a hybrid selection scheme for a crop, a broad range of genotypes could be tested in an incomplete diallel, similar to the one in [START_REF] Forst | A generalized statistical framework to assess mixing ability from incomplete mixing designs using binary or higher order variety mixtures and application to wheat[END_REF] that identified suitable material to form "pools" in cultivar mixtures. In mixed cropping, in early stages of breeding, where the size of the GMA variances of the two species are yet unknown and both species are of equal interest, an incomplete factorial with equal sizes of m and n would be advisable to subsequently design a breeding scheme based on the results. In later stages of both hybrid breeding (heterotic pools have been formed) and breeding for mixed cropping (focal species has been chosen), an incomplete factorial, e.g., in the form of design D, could be applied to both pools (hybrid breeding) or the focal species (mixed cropping). Only little literature has been published on actual experiments for breeding for mixed cropping. Some authors focus on the stepwise approach, first conducting a topcross design (similar to design C) to identify most promising genotypes for mixtures, followed by a full factorial to identify best combinations (similar to design B) for the development of two components of a mixture, such as species mixtures of maize [START_REF] Hoppe | Entwicklung von Energiemaissorten für die Mischkultur mit Stangenbohnen[END_REF] or common bean [START_REF] Starke | Selektion von Stangenbohnensorten (Phaseolus vulgaris L.) für den Mischanbau mit Mais[END_REF]. The results presented in this study suggest such stepwise experiments could have been combined to a single one by the application of an incomplete design, thus speeding up the selection process. Incomplete designs can be applied to similar problems where factorial experimental designs are used, notably hybrid breeding, animal breeding.

Pr/As Concept to Select Genotypes According to Their Species-Specific Mixing Ability

The Pr/As concept allows to cluster genotypes into groups of particular "mixing-behaviors" within positive or negative GMA. Therefore, depending on the desired proportion of fraction yield by farmers, in our example a larger ratio of pea to barley, pea genotypes can be selected from either sectors U or V of Figure 3, while pea genotypes of sector W would support a higher proportion of barley. The Pr/As concept also allows to select for GMA maximization via a regression of the Pr effects on the As effects, i.e., fit a regression to the Figure 3 dataset. For instance, a regression with a slope strictly steeper than -1 (e.g., -1.5) indicates total yield can be increased by more competitive cultivars of the focal species, thus, GMA is maximized via the selection toward higher Pr or lower As effects.

The Pr/As concept can be seen as the genetic correspondence to the replacement series as, for example, described by [START_REF] Wendling | Specific interactions leading to transgressive overyielding in cover crop mixtures[END_REF]. They compared four different crop species in pairwise combinations for their biomass yield under mixed cropping. Both the Pr/As concept and replacement series describe levels of competitiveness between two species under varying competitive conditions within the mixture, conveyed by different sowing ratios in the replacement series, and by genetic differences in competitiveness in the Pr/As concept. The low and high sowing ratios of a species in the replacement series would then correspond to low (positive As effects) or high competitive (negative As effects) genotypes, respectively. [START_REF] Wendling | Specific interactions leading to transgressive overyielding in cover crop mixtures[END_REF] observed linear relationships between mixed crop species only in two out of twelve replacement series. For all other scenarios, local maxima with transgressive overyielding were identified, i.e., mixture biomass yield exceeded the pure-stand biomass yield of each species. Due to the resemblance of the two concepts, it is quite possible, that similar local maxima for mixture yield occur in a Pr/As context. In this case, instead of a linear regression, bi-, polynomial, local, or non-parametric regressions could be applied, in order to find a target interval for As values to maximize mixture yield.

Pr and As effects are correctly estimated without being biased by different levels of correlated errors. A bivariate model is provided, able to take such a correlation into account. This model is an original analysis approach and the canonical way to treat paired variables that are presenting obvious inter-dependencies (yield of two crops cultivated on the same plot). However, this bivariate approach, with our design chosen to ensure balance and avoid confounding as much as possible, did not yield an improvement in terms of precision of estimates compared with the univariate approach. This is against our expectation and suggesting literature [START_REF] Sørensen | A comparison of bivariate and univariate QTL mapping in livestock populations[END_REF][START_REF] Meier | Joint analysis of psychiatric disorders increases accuracy of risk prediction for schizophrenia, bipolar disorder, and major depressive disorder[END_REF] but the focus of this study was on parameter estimation, whereas strengths of bivariate approaches might rather lie in other applications, such as hypothesis testing and outcome prediction. Even though the precision of estimates could not be improved with it, the multivariate approach can be used to estimate genetic correlations between traits [START_REF] Meyer | Estimating variances and covariances for multivariate animal models by restricted maximum likelihood[END_REF], which is fundamental for the use of indirect selection methods in mixed crops as suggested by [START_REF] Annicchiarico | Do we need specific breeding for legumebased mixtures?[END_REF].

Pr/As-Trait Relationships to Shape Species Mixtures

The Pr/As concept can be seen as an extension of the concept of competitive effect and response [START_REF] Goldberg | Competitive effect and response in four annual plants[END_REF]. There, a relationship or effect of trait A (e.g., early vigor), measured in species one (competitive effect) on a different trait B (e.g., yield), measured in species two (competitive response) is assumed. In the Pr/As context, however, trait A in species one (e.g., early vigor) can have an effect on a trait B that is common in both species (e.g., yield of species one and two, i.e., Pr/As effects), as visualized in Figure 4. By combining the Pr/As concept with trait measurements, the BIF of a trait can be determined. This bears the potential for further systematic investigation and categorization of trait functions in mixed cropping and community ecology, where it might serve to discover, which trait categories prevail in successful plant or other organismic communities, and shape the functioning -or non-functioning -coexistence of these. The identification of BIFs is therefore very important for breeding for crop mixture, which is not possible, if only the trait-GMA relationship (i.e., total yield) is being looked at. Nevertheless, analyzing correlations between GMA and traits can still be of interest to identify key traits that influence total mixture performance like in forage crops but not the performance of individual species.

CONCLUSION

A GMA/SMA model could be applied to compare different experimental designs for their capacity to provide meaningful information for breeders engaging in mixed cropping. Based on our findings, we recommend to use an incomplete factorial design in early stages of breeding for mixed cropping, since it allows to extend the number of tested cultivars at equal levels of experimental resources. Breeding programs can be sped up by the possibility to merge otherwise stepwise full factorial experiments into one single step. The Pr/As concept applied to an incomplete factorial design was shown to be an adequate tool to optimize mixture yield. It enables (i) to select genotypes with a suitable GMA-type thus optimizing mixture composition and (ii) to identify competitive optima for yield maximization in mixed cropping. It further allows to characterize the function of traits within species mixtures by their BIF and thus gain knowledge about their role in the biological interactions between species in a plant community. For ease of comprehension, the current study does not take genotype × environment (G×E) interactions into account. Future research should address these interactions, and the models and methodology provided here can be expanded to integrate this interaction. This study provides an integrative methodological approach for the emerging field of breeding for mixed cropping of arable crops.

FIGURE 1 |

 1 FIGURE 1 | The four experimental designs used in this study, comprising three full (f) factorials with all possible pairwise combinations present and an incomplete (i) factorial with only a subset of all possible pairwise combinations present. The designs included (A) an 8 barleys × 30 peas full factorial (8 × 30f), (B) an 8 barleys × 8 peas full factorial (8 × 8f), (C) an 2 barleys × 30 peas full factorial (2 × 30f), and (D) an 8 barleys times 30 peas incomplete factorial (8 × 30i). The last three designs consume roughly the same amount of experimental resources with 64, 60, and 60 experimental plots, respectively.

FIGURE 2 |

 2 FIGURE 2 | Decomposition of the general mixing ability (GMA) of pea G p i in its producer P p i , and associate A p i effects. Parameters apply in analogy for barley.

Figure 1

 1 Figure12 means of 1000 estimates per design and 95% confidence intervals 3 means of rho-z-rho transformed Pearson correlations, see "Materials and Methods," letters a-d

FIGURE 4 |

 4 FIGURE 4 | Schematic representation of pea genotypes with three potential relationships (positively correlated, uncorrelated, negatively correlated) of their GMA with a fictive explanatory trait and three potential underlying Pr-and As-trait relationships. Values of the explanatory trait lie on the x-axis, GMA, Pr, and As values on the y-axis. Pr/As-trait relationships reveal different biologic interaction functions (BIFs). The pattern describes a neutral (0), positive (+) or negative (-) influence on the species on which the trait was measured (left of the slash) or the species associated to this species (right of the slash).

TABLE 1 |

 1 Results on general and specific mixing ability (GMA and SMA) from simulated data for four different trial designs (A-D).

	Error	Means of variance	estimates	(truth: 5.0)	5.02 ± 0.03	5.02 ± 0.06	5.03 ± 0.00	5.02 ± 0.06	4.88 ± 0.02	4.76 ± 0.04	4.76 ± 0.05	4.80 ± 0.05
	SMA	Means of correlations	BLUPs vs. true	effects	0.77 d	0.74 c	0.70 b	0.67 a	-5	-	-	-
		Means of variance	estimates	(truth: 5.0/0.0)	5.00 ± 0.05	5.03 ± 0.10	4.85 ± 0.12	4.95 ± 0.13	0.14 ± 0.01	0.27 ± 0.02	0.35 ± 0.03	0.36 ± 0.03
	Barley GMA	Means of variance Means of correlations	estimates BLUPs vs. true	(truth: 5.0) effects	4.98 ± 0.17 0.98 c	4.98 ± 0.19 0.93 b	4.64 ± 0.43 -4	5.13 ± 0.21 0.91 a	4.94 ± 0.17 0.99 c	4.94 ± 0.18 0.98 b	4.64 ± 0.43 -4	4.89 ± 0.18 0.96 a
	Pea GMA	Means of variance Means of correlations	estimates BLUPs vs. true	(truth: 3.0) 2 effects 3	2.97 ± 0.06 0.88 c	2.90 ± 0.13 0.89 d	3.06 ± 0.12 0.68 b	3.01 ± 0.13 0.65 a	2.98 ± 0.05 0.95 c	2.90 ± 0.11 0.96 d	2.89 ± 0.07 0.84 b	2.84 ± 0.07 0.83 a
		Plots per	replicate		240	64	60	60	240	64	60	60
		Trial	design 1		A (8 × 30f)	B (8 × 8f)	C (2 × 30f)	D (8 × 30i)	A (8 × 30f)	B (8 × 8f)	C (2 × 30f)	D (8 × 30i)
		Simulation			SMA-present	simulation			SMA-absent	simulation

Frontiers in Plant Science | www.frontiersin.org

January 2021 | Volume 11 | Article 620400

ACKNOWLEDGMENTS

The authors want to thank Hans-Peter Piepho and Nathalie Moutier for valuable discussions and feedback.

FUNDING This work was a part of the project ReMIX "Redesigning European cropping systems based on species MIXtures, " funded by the EU's Horizon 2020 Research and Innovation Programme (Grant Agreement No. 727217) and the Swiss State Secretariat for Education, Research and Innovation (SERI, contract number 17.00091).

DATA AVAILABILITY STATEMENT

The R-codes used for the simulations and analyses supporting the conclusions of this article were made available by the authors under (see also Haug, 2020a,b).

AUTHOR CONTRIBUTIONS

BH wrote the manuscript. The simulation experiments were designed and its statistical framework was established by all authors. IG, JE, TMH, EF, MM, PH, and BH developed the incomplete design idea. TF, TMH, and BH wrote the R-scripts for the analyses. All authors reviewed the final manuscript.

SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2020. 620400/full#supplementary-material

Conflict of Interest:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


W Wang

A L Helbert

F Brisset

M H Mathon

T Baudin

Monte Carlo simulation of primary recrystallization and annealing twinning

Keywords: Monte Carlo simulation, primary recrystallization, twinning, nickel alloy

HAL is

Introduction

The first observation of annealing twins can be traced back to the beginning of the 20 th century [1]. With the technology and material developments, it was largely observed during the annealing of the FCC materials with low stacking fault energy. The related mechanisms have also been studied for decades. Combining the proposed mechanisms [2][3][4], the annealing twin formation was defined as a nucleation of partial Shockley loops at the {111} ledges (or emission of partial dislocations) from the migrating grain boundaries (GBs) [5][6], including the formation of fivefold twin in nanocrystalline materials [START_REF] Hofmeister | Encyclopedia of Nanoscience and Nanotechnology[END_REF]. This has been supported by molecular dynamics simulations [START_REF] Farkas | [END_REF][9]. From Gleiter's theory of [3], the twinning probability depends both on the annealing temperature and on the difference of Gibbs' free energy between the growing and the shrinking grains [3]. Pande et al. [10] found that the number of twins per grain depends on the grain size (not directly on temperature), and on the driving force. However, in these theories, the Gibbs' free energy gradient or the driving force was only related to the GB energy. They were then modified by Cahoon et al. [11] and Li et al. [12], who introduced the stored energy effect. Moreover, with the Gleiter's model [3] the annealing twin was formed as an accident when the migrating GB was parallel to a {111} plane of a growing grain. This theory was verified by a recent study based on transmission electron microscope observations (TEM) [START_REF] Wang | Etude des mécanismes de recristallisation de nouveaux alliages à base de Ni[END_REF] and shown that the twin formation depends on the direction of grain growth. All of these studies show that the annealing twinning is an important mechanism related to primary recrystallization.

As a statistical and stochastic approach, Monte Carlo (MC) simulations are largely used for primary recrystallization and grain growth. However, the twin formation has not yet been introduced in this type of simulation, especially using the orientation image microscopy (OIM™) maps as input data, which contain the microstructure, the position of each pixel and their orientation.

Before the development of field emission gun-scanning electron microscopy (FEG-SEM), it was difficult to obtain an OIM map of the deformed state due to the resolution limit. Despite this problem, Caleyo et al. [START_REF] Caleyo | [END_REF] have tried to simulate the primary recrystallization from an OIM map measured by SEM, but the non-defined pixels were attributed the orientations considering the global texture of the deformed material. Baudin et al. [15] have also run simulations using TEM data. Note that the input data could also come from finite element modelling of a deformed state [16]. After the development of FEG-SEM, the MC simulations directly from OIM map became feasible [17][18].

In the present paper, the authors will realize a microscopic simulation of the primary recrystallization through MC approach. The experimental data of a deformed material, obtained by electron backscatter diffraction (EBSD) system installed on a FEG-SEM will be used as input data.

The multiple twinning will also be introduced with considering the variation of stored energy and the direction of grain growth.

Material and initial structure

A nickel alloy after hot forging and a homogenization process was employed [19]. After these treatments, this material presents a recrystallized and homogenous microstructure. The mean grain size is about 17 µm without considering the twins.

A sample cut from this material, was deformed by cold rolling with 60% reduction. Its deformed microstructure was analyzed by an EBSD system installed on a FEG-SEM. This analysis was carried out on the transverse plane (Rolling Direction RD -Normal Direction ND).

The acquisition step was 200 nm based on a hexagonal grid. Then an EBSD map of 50 × 50 µm 2 was chosen as the initial microstructure for simulation, whose texture is in good agreement with the global texture (Figure 1). The main texture components are Copper {112} <111> and the α -fibre (from Goss {011} <100> to Brass {011} <211>). This global texture was analyzed by neutron diffraction at the Laboratoire Léon Brillouin (CEA Saclay, France) on the diffractometer 6T1.

The distribution of stored energy within this chosen microstructure was obtained, with the help of line broadening measurements of neutron diffraction peaks [20][21] and the calculation of kernel average misorientation (KAM) [START_REF] Brewer | Electron Backscatter Diffraction in Materials Science[END_REF][START_REF] Calcagnotto | [END_REF]. 1 ( , , ) ( , , )

n neutron i i i i H F H          (1) 12 ( , , ) i i i
F    is the volume fraction of every orientation calculated with the LaboTex software from the orientation distribution function (ODF). The ODF is determined by using the discrete ADC method [24]. This mean value allows us to have a global vision of the stored energy in this deformed material. Moreover, the local deformation within the microstructure was calculated by the KAM parameter obtained from EBSD data. The kernel size is chosen as the 1 st neighboring pixels. The minimum misorientation is chosen as 15° to define the GB. For every pixel, if the misorientation between it and any of its 1 st neighboring pixel is larger than this value, this neighboring pixel will be considered as belonging to another grain, and its misorientation will be removed from the KAM calculation. After this calculation, every pixel i has a KAM value i  . According to that, the pixels can be color coded to obtain the KAM map, which will show the deformation distribution, associated to geometrical necessary dislocations, within the microstructure. As shown in figure 3, the local deformation is always around the high angle grain boundaries (HAGBs) (angle > 15°) (red lines) and the distribution is almost the same in every grain. From the KAM value i  of all the pixels, the maximum max  and the minimum min  misorientation within the microstructure can be detected. Then the stored energy of every pixel is distributed accordingly to the following rules:

The

MC simulation model

In the MC simulations and for every calculation, a pixel i will be firstly chosen randomly.

The simulation will then calculate the growth probabilities of its neighboring pixels towards itself (Figure 4a) [START_REF] Landau | A guide to Monte-Carlo simulations in statistical physics[END_REF][START_REF] Rollett | Continuum scale simulation of engineering materials: Fundamentals-Microstructures-Process application[END_REF]. The recrystallization corresponds to the growth of the neighbor pixel that has the highest probability (as the first neighbor in figure 4a).

Figure 4:

(a) Typical model of MC simulation: for a chosen pixel i, the recrystallization correspond to the growth of one neighboring pixel towards itself, which has the highest probability (as the growth of the 1 st neighbor towards the pixel i); (b) our model: for a chosen pixel i, the recrystallization correspond to its growth towards one neighboring pixel which has the highest probability (as the growth of pixel i towards the 1 st neighbor).

In our simulation, the twinning will be introduced during primary recrystallization. If the typical model is used, every neighbor has the twinning probability towards the pixel i. Moreover, every neighbor has four <111> axes as the possible twinning systems. So for every calculation, 24 twinning systems (4 axes x 6 neighbors) should be considered. In our model, the recrystallization corresponds to the growth of one chosen pixel i towards one of its neighbors with the highest possibility (as the first neighbor, figure 4b). Therefore, only the 4 axes of the chosen pixel i will be considered as the twinning systems.

In addition, the MC simulation was used only for the primary recrystallization and associated twinning during annealing at 900°C. Indeed, the experimental results were shown in a previous paper [19], and the recovery is considered negligible during annealing. Therefore, the recovery is not introduced in the present study. More, the grain growth after complete primary recrystallization is not simulated here.

Primary recrystallization model

The primary recrystallization is simulated in two parts: nucleation and growth of formed nuclei.

Nucleation

The recovery is an important process for the subsequent nucleation, as well as the time for nucleation and the distribution of nuclei. However, in our simulation, there is no recovery process.

Instead of it, the critical number of nuclei nuclei CN is used as a condition of nucleation. As a new recrystallized grain directly provides from the growth of a nucleus, the number of nuclei can be estimated experimentally from the number of grains in the sample after complete recrystallization.

The critical nucleus size is about 200 nm, observed experimentally by TEM [START_REF] Wang | Etude des mécanismes de recristallisation de nouveaux alliages à base de Ni[END_REF]. As this value corresponds to the acquisition step size of an EBSD analysis, every pixel in the initial microstructure can be considered as a possible nucleus. To be a real nucleus, there are two critical conditions on this pixel:

 It depends on the probability of its growth towards one of its neighboring pixels (explained in the next section);

 The number of existing nuclei is still lower than the critical value nuclei CN .

Growth of the nuclei

During the simulation, the growth of a nucleus corresponds to its growth towards the deformed area. The method described in [START_REF] Volovitch | [END_REF] was used with some modifications. For a chosen pixel i (belonging to a nucleus), the driving force for its growth depends on the stored energy

gradient ij H  , the GB energy () ij
 and the mobility

() ij M  ( ij
 is the misorientation between the pixel i and its j th neighbor). The probability ij P for the growth of pixel i towards its j th neighbor is given by the following equations 3-4:

( ) [ ( ) ] ij ij ij ij P M C H       if 0 ij E  (3) ( ) [ ( ) ] exp( / ) ij ij ij ij ij P M C H E kT        if 0 ij E  (4)
C is a constant to show the ratio between the influence of the stored energy and that of the GB energy during primary recrystallization. Its value is in the range of 0.1 and 10 [START_REF] Volovitch | [END_REF]. In the present paper, both influences were supposed to have the same level (C=1).

ij E 

is the energy variation after the growth of pixel i towards its j th neighbor (or reorientation of its j th neighbor to have the same orientation as i). It can be presented by the variation of GB energy ij EG  and that of the stored energy ij ES  (Equation 5).

ij ij ij E EG ES      (5) ij EG 
is related to the misorientation variation between the pixels. After the growth, the misorientations between reoriented pixel (as 1 st neighbor in figure 4b) and its own neighbors (as A, B, C, i, 2 nd and 6 th neighbors in figure 4b) will be changed. The GB energy () jk  will be hence changed (Equation 6). For this reason, during the calculation of ij EG 

, some second order neighbors of pixel i should be considered (as A, B, C in figure 4b). and its own neighbors, respectively before and after the growth of pixel i towards j th neighbor. S is the surface of the interface between two pixels. It is constant and will be explained later.

6 6 1 1 [ ( ) ( )] final initial ij jk jk k k EG S           (6)

ij ES 

corresponds to the stored energy variation of the reoriented pixel (Equation 7). For the pixel i which belongs to an existing nucleus, its stored energy is already zero, and the stored energy variation equals to that of the reoriented pixel j.

[ ]

final initial ij j j EG H H V    ( 7
)
initial j H and final j H are the stored energies (J.mol -1 is modified to be J.m -3 ) of reoriented pixel before and after the growth of pixel i. V is the pixel volume which is constant.

The pixels are supposed to be as three dimensional ones with the acquisition step size considered as the pixel height (Figure 5). Then, the common surface S (red plane) between two pixels and the pixel volume V can be estimated.

ij     ;  If 55 1 ij    , 3 0 ( ) 1 exp[( ) ] ij ij M       with 0 15    (8)
and () ij  is calculated by the Read and Shockley equation [START_REF] Read | Dislocations in crystals[END_REF]:

max 00 ( ) [1 log( )] ij ij ij       (9)
Moreover, the GB energy and mobility for the coherent twin boundaries Σ3coh are almost zero. On the contrary, those of incoherent ones Σ3incoh are not null. However, the determination of the coherence of Σ3 GBs (including Σ3coh and Σ3incoh) is feasible but very heavy [START_REF] Wright | [END_REF].

Therefore, these two parameters are set to zero for all the Σ3 GBs. The Σ3 GBs can be easily determined according to the Brandon's definition (60° about a common <111> crystal axis with 8.7° as the maximum deviation [31]).

After the calculation of the growth probability ij P , the maximum probability will be chosen as the growth direction of the pixel i. If several neighbors have simultaneously the maximum probability, the pixel i will grow along one direction chosen randomly among them.

Twinning

Twinning with formation of Σ3coh boundary

From Gleiter's theory [3] and the TEM observation [START_REF] Wang | Etude des mécanismes de recristallisation de nouveaux alliages à base de Ni[END_REF], the active twinning system seems to depend on the direction of grain growth. Therefore, in this simulation, the growth direction of a growing grain will be firstly detected according to the process described in section 3.1 (blue arrow towards 6 th neighbors in figure 6a). Then, the directions of its 4 <111> axes will be also analyzed (green lines in figure 6a). If there is at least one <111> axis that approaches the growth direction (dotted green line), it is possible to have twinning from this growing grain. If more than one <111> axis is aligned with this direction, the closest one will be taken as the twinning system. In addition, during the annealing twin formation, a new Σ3coh boundary will be formed. The GB Σ3incoh are formed as the indirect outcomes after the formation of coherent ones Σ3coh [START_REF] Wang | Etude des mécanismes de recristallisation de nouveaux alliages à base de Ni[END_REF].

There are two conditions for a boundary to be considered as a Σ3coh boundary between two grains:

 The misorientation of this interface should correspond to a 60° rotation around a common <111> axis (as the definition of Σ3) and  The plane of the Σ3coh boundary should be parallel to a common {111} plane of both grains, or perpendicular to the common <111> axis of rotation.

To state on the coherence of Σ3 boundary, both previous conditions must be checked. The first condition can be verified by the twin orientation calculation (as in figure 6b, 6 th neighbors is reoriented with the twin orientation of pixel i) [32]. For the second condition, the orientation of the boundary plane should be considered and compared to a {111} plane. Nevertheless, this is not easy to calculate the boundary plane (or plane trace) from an EBSD map. As shown in figure 6b, the new formed Σ3coh boundary (red line) is not perpendicular to the common <111> axis, and it deviates thus from the ideal Σ3coh direction (orange line). This deviation is induced by the hexagonal pixel shape.

For the EBSD map which contains hexagonal pixels, the grain boundaries can be defined by OIM in two ways: either a serie of segments between the adjacent pixels (red line in figure 7) or a straight line between two triple points (blue line in figure 7). In this simulation, the new formed GB is as short as one segment between two pixels (Figure 8a). Both GB definitions of OIM are the same, and this short GB has only three possible directions (as a, b and c shown in figure 8a). But, for the real situation, the ideal Σ3coh boundary could be oriented towards any direction. Consequently, it is almost impossible to simulate the ideal Σ3coh boundary with limited GB directions at the first stage of twin formation. For example, if the GB corresponds to three segments (Figure 8b), both definitions become different. The number of possible GB directions increases. By analogy, the more the GB is long, the more possible directions have the GBs [START_REF] Wright | [END_REF]. Then, it becomes easier to have the GBs parallel to the <111> axes.

Therefore, in our simulation, the twinning will be firstly realized according to the <111> axis of parent grain that is close to its growth direction, with the formation of a Σ3coh boundary (black line in figure 9b). During the growth of the parent grain and its twin, the length of the formed Σ3coh will increase. At the same time, the orientation of the Σ3coh will be adjusted by comparison to the ideal one (Figure 9c-d). In this figure, the Σ3coh boundary is defined by the second method presented in figure 7.

Multiple twinning

In fact, during the experimental annealing, the twin formation corresponds to multiple twinning. It is revealed by the presence of several Σ3coh boundaries or twin parts in the same recrystallized grain. These twin boundaries are either parallel to (retwinning [32]) or intersecting with each other.

The stored energy effect on twin formation has been already observed in the same alloy [19].

It was found that the twinning depends on the stored energy E (or strain amount) and the grain size D (Equation 10).

0 ( )( ) N A E B D D    (10)
A and B are constants. N is the Σ3 boundary (including Σ3coh and Σ3incoh) number per grain. 0 D is the critical grain size for the twin formation.

For multiple twinning, a simplified equation 10 was employed (Equation 11).

' ( ) N A D D  (11)

Indeed, since the simulation is realized for only one deformation rate, (AE+B) in replaced by the constant A' , and N is considered as ideal Σ3coh boundary number or twinning times of a grain.

Then the twin formation was only controlled by the grain size. For a grain of size D, if the existing Σ3coh boundary number is smaller than the experimental value N, the new twinning is authorized during its growth.

Results and discussion

Recrystallization

This simulation is used to simulate the recrystallization and twinning of a deformed alloy during an annealing at 900°C. Experimentally, the recrystallization is very fast, and the recovery is negligible [19].

Nucleation

The figure 10

Kinetic of recrystallization

The recrystallized fraction during simulation is plotted in figure 11, as a function of the MC step (on log scale). From the Johnson-Mehl-Avrami-Kolmogorov (JMAK) equation, the Avrami exponent n is 2. This simulated value is close to the theoretical one which should be in the range between 2 and 3 for a two dimensional growth [START_REF] Humphreys | Recrystallization and related annealing phenomena[END_REF]. This theoretical value range is obtained on the basis of a homogenous distribution of nuclei during recrystallization. The experimental value is always smaller than 2 [START_REF] Humphreys | Recrystallization and related annealing phenomena[END_REF]. This is due to the heterogeneous distribution of nuclei, as effect of shear band which is a favor site for nucleation. For our simulation, the nuclei distribution is homogeneous (Figure 10). During experimental annealing, the distribution of nuclei is heterogeneous with formation of nucleus clusters. Their growth will be blocked by each other during the recrystallization because of the null stored energy gradient. This growth restraint delays the recrystallization kinetic. Therefore, the bigger value of n compared with experimental one is reasonable. Moreover, from JMAK equation, the pre-exponential factor k is 0.0049 which depends on the temperature [START_REF] Avrami | [END_REF]. However, it is difficult to compare this simulated parameter with experimental one. Because the recrystallization during simulation is in function of MC steps, but that during experiment is related to time.

Microstructure after complete recrystallization

The recrystallized microstructures after simulation and an annealing of 4 min at 900°C (experiment) are shown in figure 12. Within these microstructures, the colors are used only to separate the grain without considering the twin boundaries. One can find that the simulated microstructure is similar to the experimental one. However, the grains within experimental microstructure seem to be bigger than simulated ones. Therefore, more statistical analyses as the mean grain size and the grain size distribution were carried out for the experimental microstructure. These analyses are done from an EBSD map of 150 x 350µm 2 for a statistical result. Nevertheless, the same analyses were also done for the simulation, from the total simulated area of 50 x 50 µm 2 . The mean grain sizes are 4.7 and 4.1 µm for simulation and experiment, respectively. Then, the mean grain size is almost identical for both microstructures. The grain size distributions of both microstructures are shown in figure 13. From the distribution curves, there is a difference between the peak values of both distributions. The experimental microstructure presents more small grains and more coarse grains (black arrows) than simulated one. One can assume that after complete recrystallization, some grains grow by consuming some others that become small. This grain growth is not introduced in the simulation, and then the recrystallized grains cannot thus grow at the expense of others.

Twinning

Morphology of twin boundaries

Four morphologies of annealing twins and thus the twin boundaries (Figure 14a) have been presented by Mahajan et al. [5]. A is a twinned area at the corner of the parent grain with one Σ3coh boundary (red line). B is a twin traversing the grain with two parallel Σ3coh boundaries. C is a twin which terminates in the grain. It is surrounded by two Σ3coh boundaries and a Σ3incoh (green line). As the twin formation is three dimensional, the situation D which locates inside of the parent grain is also possible. It corresponds to the trace of a twin. Besides these situations above, two other morphologies were also observed by MET studies of the same alloy (Figure 14b andc) [START_REF] Wang | Etude des mécanismes de recristallisation de nouveaux alliages à base de Ni[END_REF]. E is a twinned part that traverses the grain. However, the interface between it and the parent grain correspond to two Σ3coh linked by a Σ3incoh. It was also found that the parent grain can twin towards different directions with formation of two twinned parts as F and G. When they meet together during their growth, a Σ9 (blue line) is formed.

For this situation, the Σ3coh boundaries are intersecting and so are not parallel. From this figure, all the experimental morphologies shown in figure 14 are reproduced by simulation. The figure 16a shows the situation B. The figure 16b presents the morphology F-G with two Σ3coh and one Σ9 from a triple point (blue arrow). The figure 16c shows the simulations of twins A and C (orange arrow reveal the Σ3incoh boundary). Finally, the situations E and D are also obtained (Figure 16d). As only the Σ3coh boundary formation is introduced in this simulation, the presence of the Σ3incoh and Σ9 boundaries is an indirect outcome. Indeed, the Σ3incoh and Σ9 boundary formation corresponds to the association of two Σ3coh boundaries. This phenomenon is largely studied during grain boundary engineering (GBE) [35]. From this simulation, the formation of situation E (Figure 14b) can be explained by figure 17a. During the growth of parent grain 1, one side (on the right) forms a twinned part E and creates a Σ3coh boundary (green line). On the contrary, the other side (on the left) of parent grain 1 forms a twin later and at another location. When the two twinned parts E meet together, the Σ3incoh boundary is formed. Therefore, the situation E can be considered as the result of two Σ3coh formed at different moment or different location.

The situation C with presence of a Σ3incoh boundary is similar to the situation E. The schema of this case is shown in figure 17b. During the growth of parent grain 1, the twinned part C is formed from the right side as situation E with creating a Σ3coh boundary. On the contrary, instead of twin formation, the other side (on the left) of parent grain 1 continues to grow. At the same time, the twinned part C grows also. During its growth, the area C twins and develops another parallel Σ3coh boundary. The Σ3incoh boundary is also formed automatically as a geometrical product. This retwinning process allows the formation of an area (parent grain 2) with the same orientation as parent grain 1. During the growth of both parent grain parts, they meet each other and a closed "circle" grain is formed. Therefore, the situation C is the result of twice twinning at one side of parent grain and the growth of parent grain surrounding the twinned part C.

From situations E and C, one can find, the Σ3incoh boundary was well simulated as an indirect product which is related to formation of Σ3coh ones. The situation D in figure 14a is considered as the trace of a twin part, because the twin formation is three dimensional. However, our simulation is two dimensional, so the formation of this situation (twin part T' indicated by black arrow) as shown in figure 16d is not physical. In fact, the area T' is formed fist by twinning from T (straight boundary). Then the growth of T surrounds it with formation of a closed area T'. This simulated twin T' is surrounded by a straight Σ3coh boundary and a semi-circle Σ3incoh boundary. Let us note, that as these two boundaries form a closed "circle", they cannot be separated by the OIM analysis and the coherence of the Σ3 cannot thus be detected (Figure 16d).

Twinning frequency

The twinning frequency is studied by analysis of Σ3 (including Σ3coh and Σ3incoh) fraction and its density. The definition of these two parameters has been already presented in [19]. The corresponding values obtained from simulation and experiment after 4 min annealing at 900°C are shown in table 1. This is because the twinning during simulation is controlled by two conditions:

 the evolution of twin formation (Equation 11) which is obtained through experiments;

 the presence of at least one <111> crystal axis of parent grain close to its growth direction, itself chosen by probability calculation.

For the experimental situation, as shown in figure 18a, the grain 1 can grow towards grain 2 simultaneously in all directions with the same probability (red arrows). Then it is easy to have one <111> axis (green lines in pixel i) that is parallel to the growth direction for twin formation. On the contrary, during the MC simulation, only one growth direction is chosen. It is thus possible to have none of corresponding <111> axis (Figure 18b). Therefore, the used model decreases the twin frequency during simulation.

Conclusions

The MC simulation of primary recrystallization at 900°C was realized for a deformed nickel alloy. Without recovery process, the simulated nucleation is homogeneous. This favors an Avrami exponent n close to the ideal value, but deviated from the experimental one. Moreover, the recrystallization kinetic is in agreement with experimental results. After complete recrystallization, because the grain growth is not introduced in the simulation, a difference of grain size distribution between simulation and experiment is observed. However, the mean grain size are similar.

In addition, the annealing twinning was also simulated at a microscopic scale. Different morphologies observed during experiments were simulated. Therefore, the annealing twin process can be considered as a mechanism of primary recrystallization. The dependence of grain growth direction on twin formation reveals the possibility of parallel Σ3coh boundaries formation.

Moreover, it showed the Σ3incoh and Σ9 boundaries are the indirect outcome after Σ3coh formation. A smaller twin frequency is obtained during simulation compared with experiment because at this time, only one grain growth direction is taken into account in the simulation.

Figure 1 :

 1 Figure 1: (a) Bulk texture (neutron diffraction) of the deformed sample after 60% cold rolling, and (b) local texture (OIM) of the initial microstructure for the simulation.

Figure 2 :

 2 Figure 2: Stored energy of the main components after cold rolling with 60% reduction.

Figure 3 :

 3 Figure 3: (a) EBSD map (Image Quality -IQ) of the initial microstructure and (b) KAM map of the same area. The red lines show the HAGBs.

  energies (J.m -2 ) between reoriented pixel (j th neighbor)

Figure 5 :

 5 Figure 5: Supposed three dimensional pixels with acquisition step size a considered as the pixel height. The red plane is the common surface between two adjacent pixels.

Figure 6 :

 6 Figure 6: Twinning model: (a) comparison of <111> crystal axes (green lines) with chosen growth direction (blue arrow) and (b) twinning with the axis the most close to growth direction (dotted green line). The orange line indicates the theoretical Σ3coh orientation which is perpendicular to the twinning axis. The red line shows the simulated Σ3coh after twinning.

Figure 7 :

 7 Figure 7: GB definitions by OIM. The red line shows the definition with a series of pixel segments. The blue line presents the GB as a straight line between two triple points.

Figure 8 :

 8 Figure 8: Short grain boundaries (one segment) and long grain boundaries (three segments) from a triple point A. the red lines and blue lines show respectively two GB definitions.

Figure 9 :

 9 Figure 9: An example for (a) a recrystallized grain (red pixel) surrounded by the deformed area (purple); (b) twinning with formation of a short Σ3coh boundary; (c-d) growth of parent grain and its twin (blue pixel) with adjustment of Σ3coh boundary compared with the ideal one (orange line). In this figure, the Σ3coh boundary is defined by the second method presented in figure 7.

  presents the nucleation after one step of MC simulation. It shows the nuclei are formed close to the area with strong local deformation (black areas). The distribution of these nuclei is homogenous within the deformed microstructure. This homogenous distribution of nucleus is due to the lack of recovery process during simulation which controls the location and moment of nucleation.

Figure 10 :

 10 Figure 10: IQ map for the distribution of nuclei after one step of MC simulation (nuclei are shown by red points).

Figure 11 :

 11 Figure 11: Kinetic of recrystallization obtained by MC simulation.

Figure 12 :

 12 Figure 12: Microstructure after (a) simulation and (b) annealing of 4 min at 900°C. The colors are used only to distinguish the grains without considering the twin boundaries.

Figure 13 :

 13 Figure 13: Grain size distribution after simulation and annealing of 4 min at 900°C (experiment).

Figure 14 :

 14 Figure 14: (a) Four different morphologies of twin parts presented by Mahajan et al. [5], (b) and (c) two other morphologies found during MET study of the same alloy. The red, green and blue lines show respectively, the Σ3coh, Σ3incoh and Σ9 boundaries.

Figure 15 :

 15 Figure 15: Distribution of GBs within the recrystallized microstructure. The black, red and green lines show the HAGBs, Σ3 and Σ3coh boundaries (detected by OIM TM ).

Figure 16 :

 16 Figure 16: Detailed morphologies of twinned parts obtained by simulation. The black, red and green lines show the same GBs as figure 15. The blue and orange arrows show Σ9 and Σ3incoh, respectively.

Figure 17 :

 17 Figure 17: Schemas for the formation of situations (a) E and (b) C, shown in figure 14.

Figure 18 :

 18 Figure 18: (a) Real growth simultaneously in several directions with the same probability and (b) MC simulation with only one chosen direction. The blue lines show the HAGBs. The green lines

Table 1 :

 1 Σ3 fraction and its density obtained by simulation and experiment Both parameters reveal the twinning frequency is lower during simulation than experiment.

		Σ3 fraction	Σ3 density (µm -1 )
	Simulation	30.7%	0.37
	Experiment	51.5%	0.54

Acknowledge

The authors would like to thank Denis SOLAS, Isabelle DROUELLE and Thierry AUGER for fruitful discussions, and thank Aubert & Duval for providing materials.