domain: stat sdv

text:
Michel Turbet Delof
email: michel.turbet-delof@inrae.fr

Pierre Rivière

Julie C Dawson

Arnaud Gauffreteau

Isabelle Goldringer

Gaëlle Van Frank

Olivier David

Bayesian joint-regression analysis of

Keywords: decentralized participatory plant breeding, bread wheat, G × E interaction, hierarchical model, Finlay-Wilkinson model, Student's t distribution, varietal stability

come

Introduction

Developing new varieties adapted to Organic Agriculture (OA), agroecological and low input systems is a major concern to achieve improvements in agricultural sustainability [START_REF] Wolfe | Developments in Breeding Cereals for Organic Agriculture[END_REF]. In OA, the use of synthetic inputs (nitrogen, phytochemicals) is not allowed, therefore, cropping environments are not standardized by inputs and varieties grow in more diverse conditions from farm to farm [START_REF] Dawson | Decentralized Selection and Participatory Approaches in Plant Breeding for Low-Input Systems[END_REF].

These environments are dependent on pedoclimatic conditions, yearly weather, farmers' management practices and interactions between these factors [START_REF] Desclaux | Changes in the Concept of Genotype × Environment Interactions to Fit Agriculture Diversification and Decentralized Participatory Plant Breeding: Pluridisciplinary Point of View[END_REF].

In order to develop varieties adapted to such a diversity of environments two strategies can be used: (i) centralized and indirect selection, or (ii) decentralized and direct selection. The key difference between these approaches is the way they take genotype-by-environment (G × E) interactions into account. These interactions are considered by plant breeders as the main factor limiting the efficiency of the response to selection in breeding programs [START_REF] Ceccarelli | Farmer Participation in Barley Breeding in Syria, Morocco and Tunisia[END_REF]. In centralized and indirect selection, breeding lines are evaluated and selected at a few research stations assumed to represent the target environments. This is efficient if there is a high additive genetic correlation between the trait measured on the station and the same trait measured in the target environment, and if the narrow sense heritability is high in the selection environment [START_REF] Falconer | Introduction to Quantitative Genetics[END_REF].

Decentralized selection can take account of G×E interactions that are important in OA [START_REF] Dawson | Decentralized Selection and Participatory Approaches in Plant Breeding for Low-Input Systems[END_REF][START_REF] Murphy | Evidence of Varietal Adaptation to Organic Farming Systems[END_REF]. In this approach, the selection and evaluation environments are very close to the target environments (the production environments of farms). Selection then maximizes the use of the reproducible part of G × E interactions to select for specific adaptations [START_REF] Annicchiarico | Response of Common Wheat Varieties to Organic and Conventional Production Systems across Italian Locations, and Implications for Selection[END_REF]. This method is close to direct selection and has been shown to be effective [START_REF] Annicchiarico | Response of Common Wheat Varieties to Organic and Conventional Production Systems across Italian Locations, and Implications for Selection[END_REF][START_REF] Ceccarelli | Farmer Participation in Barley Breeding in Syria, Morocco and Tunisia[END_REF][START_REF] Murphy | Evidence of Varietal Adaptation to Organic Farming Systems[END_REF][START_REF] Smith | Participatory Plant Breeding with Maize in Mexico and Honduras[END_REF][START_REF] Virk | Increasing the Client Orientation of Maize Breeding Using Farmer Participation in Eastern India[END_REF].

Many PPB programs have been carried out over the last 20 years targeting low-input farming systems in the Global South and also OA and agroecological systems in Europe and North America (Ceccarelli and Grando, 2020). A few programs tested different experimental designs and specific statistical methods to analyze data taking GxE into account [START_REF] Mohammadi | Identifying Superior Rainfed Barley Genotypes in Farmers' Fields Using Participatory Varietal Selection[END_REF][START_REF] Snapp | Farmer Preferences and Legume Intensification for Low Nutrient Environments[END_REF]. However, few had an extensive dataset with a large number of farms, years and genotypes to allow investigation of the best relevant methods in detail. One program with such data is a wheat PPB program that started in France in 2005, as a collaboration between INRAE GQE-Le Moulon and the Farmers' Seed Network (Réseau Semences Paysannes, RSP). This PPB program had three objectives: (i) develop population-varieties adapted to farmers' practices and needs (organic management, artisanal bread quality ...) using a participatory approach, (ii) develop strategies for preserving genetic diversity through on-farm dynamic management and breeding, and (iii) learn from and improve farmers' individual and collective breeding methods and diffuse successful methods broadly.

In this program, a large number of populations was evaluated over a large network of farms in the RSP [START_REF] Dawson | Collaborative Plant Breeding for Organic Agricultural Systems in Developed Countries[END_REF][START_REF] Goldringer | Agronomic Evaluation of Bread Wheat Varieties from Participatory Breeding: A Combination of Performance and Robustness[END_REF]Rivière et al., 2015a,b;[START_REF] Van Frank | Genetic Diversity and Stability of Performance of Wheat Population Varieties Developed by Participatory Breeding[END_REF]. Because of the extensive trial network, the assessment of population performance within an environment could potentially be improved by taking account of the average performance of populations over the network and the stability of their performance, in particular temporal stability, as it determines agronomic and economic risks. Two types of stability assessment have been developed: static and dynamic stability (Becker and Leon, 1988;[START_REF] Lin | Stability Analysis: Where Do We Stand?[END_REF].

The farmers involved in the program chose which populations to evaluate on their farm, based on prior information about the parents, and characteristics of interest. As very few populations were present in all the trials, the resulting series of trials was very unbalanced, so that the estimation of population average performances and stabilities was difficult. Joint regression is a robust method for estimating genotype main effects and stability with incomplete datasets [START_REF] Finlay | The Analysis of Adaptation in a Plant-Breeding Programme[END_REF][START_REF] Pereira | Robustness of Joint Regression Analysis[END_REF]. It is based on the Finlay-Wilkinson (FW) model, which is parsimonious since the interaction effect between a genotype and an environment is modelled as the product of a genotype stability parameter, called sensitivity, and the environment main effect. Various Finlay-Wilkinson models have been used in a frequentist framework, which include fixed-effect models and models with random environmental effects [START_REF] Nabugoomu | Analysis of Series of Trials Where Varieties Differ in Sensitivity to Locations[END_REF][START_REF] Ng | Joint-Regression Analysis for Incomplete Two-way Tables[END_REF][START_REF] Patterson | Statutory and Recommended List Trials of Crop Varieties in the United Kingdom[END_REF]. In the latter models, environmental effects are assumed to come from a common distribution, thereby leading to shrunk estimates. FW models in which genotype main effects, environment main effects and genotype sensitivities are all random effects have recently been developed that can handle unbalanced data and take into account the similarity between some genotypes and the similarity between some environments. These have been implemented in a Bayesian framework and when they include random effects, these are called hierarchical models [START_REF] Carlin | Bayesian Methods for Data Analysis[END_REF][START_REF] Robert | The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation[END_REF]. Thus far, these models have been used to analyze slightly unbalanced trials [START_REF] Lian | FW: An R Package for Finlay-Wilkinson Regression That Incorporates Genomic/Pedigree Information and Covariance Structures Between Environments[END_REF]. Hierarchical joint regression has also been used to analyze very unbalanced simulated data [START_REF] Van Frank | Influence of Experimental Design on Decentralized, on-Farm Evaluation of Populations: A Simulation Study[END_REF]. This simulation study has shown that genotypes should be tested in sufficiently many trials in order to estimate their main effects and sensitivities reliably. However, this method had not been used to analyze real and very unbalanced trials. Thus, it was not clear if it could cope with the actual levels of unbalanced data seen in the French PPB on-farm trials and what insight it could give into the behavior of genotypes across environments.

Extreme data is an important issue in data analysis. In multi-environment trials (MET), they may come from either (1) errors between scoring and data formatting (measurement error, wrong labelling, etc.), or (2) environmental heterogeneity in the trial (weed infestation, soil fertility, etc.), or (3) particular environmental conditions that fall outside the normal range of environments under study (poor emergence, extreme weather, strong pest/disease pressure, etc.). In our PPB program, as cultivation environments are less controlled, extreme observations (types 2 and 3) could be more frequent than expected. This could bias estimates based on the normal distribution. Extreme observations could be removed from the dataset to solve this problem, but it is difficult to decide which observations to remove. If too many extreme observations are removed, then the variability of the data may be underestimated and the precision of the statistical analysis overestimated.

Alternatively, statistical methods that are robust to extreme observations may be used [START_REF] Hampel | Robust Statistics: The Approach Based on Influence Functions[END_REF][START_REF] Huber | Robust Statistics[END_REF]. Various robust methods have been developed in a frequentist or a Bayesian framework, in particular methods consisting in replacing the normal distribution by a Student's t distribution in statistical models. This distribution is more robust to extreme observations than the normal distribution, because it has heavier tails [START_REF] Carlin | Inference for Nonconjugate Bayesian Models Using the Gibbs Sampler[END_REF][START_REF] Choy | Scale Mixtures Distributions in Statistical Modelling[END_REF][START_REF] Lange | Robust Statistical Modeling Using the t Distribution[END_REF][START_REF] Rosa | Robust Linear Mixed Models with Normal/Independent Distributions and Bayesian MCMC Implementation[END_REF]. It has been used to handle the extreme observations of a single trial in a Bayesian framework (Besag and Higdon, 1999;[START_REF] Cao | Bayesian Inference of Spatially Correlated Random Parameters for On-Farm Experiment[END_REF][START_REF] Gianola | Prediction of complex traits: robust alternatives to best linear unbiased prediction[END_REF]. However, to our knowledge, it has not been used to analyze an unbalanced network of trials.

Recently, participatory variety trials using crowdsourcing has been used in several countries with great success [START_REF] Van Etten | Crop variety management for climate adaptation supported by citizen science[END_REF]. These methods typically use an experimental design called a triadic comparison of technologies (tricot), followed by an analysis of variety ranks (Beza et al., 2017). In the tricot design, large numbers of farmers each compare three variety subsets from the complete set of entries, and provide direct comparison rankings among them for a few traits (i.e. best/middle/worst). By using ranking methods and structuring the entry distribution as an incomplete block design, this allows for comparisons of larger numbers of varieties without overburdening individual farmers. These design options enhance breeders' ability to engage farmers in trialing experimental lines, since on-farm trials are often limited by space and farmers' time. Trialing a few experimental lines, including a check line or variety that is replicated across sites is more realistic for farmers than implementing a fully replicated design. Triadic methods are very useful in many situations, but they are not applicable to more mature networks of farmer breeders such as we have in our PPB program for wheat. The farmers in this network have selected populations over time according to their own rationale, and the populations are not randomly assigned to farmers. In addition, farmers test different numbers of populations, with some only trialing a few and others trialing several dozen. Farmers also want access to quantitative data rather than simple ranks, and so a non-parametric ranking of varieties with no assumptions about distribution will not produce a satisfactory analysis for this purpose. This article was aimed at improving the assessment of the population-varieties of our program by using the information at the level of the network. As our dataset was very unbalanced and could include extreme observations, we compared several Finlay-Wilkinson models, in particular hierarchical models and models based on the t distribution. These models were developed in a Bayesian framework, since this framework is rigorous and since it facilitates the implementation of complex models [START_REF] Carlin | Bayesian Methods for Data Analysis[END_REF][START_REF] Robert | The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation[END_REF]. Finally, the best Finlay-Wilkinson model we obtained was used to analyze our data and characterize the behaviour of our population-varieties across environments. In our study, we will call a germplasm any biological entity whose individuals are derived from the same breeding process, including varieties registered in the official catalog, landraces, historic varieties, mixtures or populations stemming from crosses. An environment is the combination of a farm and a year.

Material and methods

Germplasm

We studied 206 germplasm covering different "germplasm types": 98 "cross" germplasm resulting from crosses made either on the farm or at the research station (Rivière et al., 2015b), 50 "landraces", i.e. population varieties grown before 1850, 30 "historic varieties", developed by professional breeding before 1950, 17 "mixtures", which were generally complex, with numerous genotypes from potentially all the other germplasm types. In addition, 11 "registered varieties" after 1950 and widely used in organic farming were included:

Maitre Pierre (1954), Poncheau (1956), Renan (1990), Ataro (2004), Pollux (2004), Rubisco (2012), Hendrix (2012), Kampmann selected in Renan, and Hermes (1982), Alauda (2004) and Goldritter (2013), all three selected in Probus (1957).

Experimental designs and data

Experimental designs

Data were collected between 2008 and 2019. The wheat PPB program followed numerous experimental designs due to the different constraints of farmers, collectives and researchers. The designs have been grouped into 5 classes (Tab. 2).

The "regional farm" and "satellite farm" designs were co-designed to be adapted to the farmer-breeders' constraints and to be used in their agricultural routine. In these designs, the germplasm common to all farms (control germplasm) were collectively chosen by farmers and researchers, while each farmer individually chose the additional germplasm to be cultivated in his farm. At the beginning the control was a selection in a landrace, and after 2014 it was a germplasm stemming from a cross. Both complete-block and incompleteblock designs were used to address specific research questions such as the study of the evolution of traits (Rivière et al., 2015b), local adaptation (van Frank et al., 2020) or the evaluation of agronomic performance [START_REF] Goldringer | Agronomic Evaluation of Bread Wheat Varieties from Participatory Breeding: A Combination of Performance and Robustness[END_REF]. Some unreplicated trials corresponded to trials with replications but for which measurements could not be performed in some replications. combination of a year and a location.

Nb of blocks

Data collected

Four traits were studied, plant height (60% of the data was the average height of 25 individuals and 40% was the overall height of the microplot, mm), spike weight (mean of 25 individual measures, g), protein content of the grain (on the microplot, measured with NIRS technology at INRAE Clermont-Ferrand France, %) and thousand kernel weight (TKW, measured on the microplot, g). These four traits were among those collectively chosen by farmers and researchers to be measured during the PPB program (Table 3). Plant height was measured in the field, while the other traits were measured after harvest at the research station on samples of spikes sent by farmers. The data analyzed were the adjusted means for block effects if these effects were significant, and the empirical means if otherwise. Obvious outliers were excluded.

van [START_REF] Van Frank | Influence of Experimental Design on Decentralized, on-Farm Evaluation of Populations: A Simulation Study[END_REF] analyzed the sensitivity of the hierarchical FW model to different MET set-ups with simulated data. They found that, in contrast to the environmental effects, the germplasm effects and FW coefficients were difficult to estimate. This is why they recommended that a large number of environments be used and that the germplasm be repeated sufficiently. We have therefore made a selection of the data and kept the environments with at least five germplasm and the germplasm that were present in at least four environments. Thus, the data analyzed comprised 70 to 76% of the initial data, depending on the trait.

The multi-environment data were very unbalanced, with most of the germplasm occurring in a limited number of environments (the median number of replicates across environments was seven, and about 20% of the germplasm were replicated in four environments only). For each trait, the number of observations was between 1300 and 2000 and the measures were spread over more than nine years (Tab. 3).

Models

The phenotypic value Y ij ∈ R for a given trait Y , germplasm i and environment j was assumed to be equal to

Y ij = µ ij + ε ij ,
where (i, j) ∈ C, C was the set of the germplasm x environment combinations occurring in the data set, µ ij ∈ R was an expectation term, and ε ij ∈ R was a residual term. Five models were developed, which modelled the expectation term, the residual term and the prior distribution differently (Tab. 4).

Expectation term

In models ADHs and ADHn, the expectation term was modelled as additive effects of both the germplasm and the environment without interaction:

µ ij = α i + θ j ,
where α i ∈ R was the main effect of germplasm i, and θ j ∈ R was the main effect of environment j. Models FWHs, FWs and FWHn modelled genotype-environment interaction using the Finlay-Wilkinson, also called joint-regression, model [START_REF] Finlay | The Analysis of Adaptation in a Plant-Breeding Programme[END_REF]. In these models, the expectation term was assumed to be equal to

µ ij = α i + θ j + η i θ j ,
where η i ∈ R was the sensitivity of germplasm i to environments (linear regression coefficient, [START_REF] Perkins | Environmental and Genotype-Environmental Components of Variability[END_REF]. [START_REF] Finlay | The Analysis of Adaptation in a Plant-Breeding Programme[END_REF] defined their coefficient as b i = 1 + η i . As the average sensitivity is equal to 0, a germplasm with η i > 0 (resp. η i < 0) is more (resp. less) sensitive to environments than a germplasm with the average sensitivity [START_REF] Nabugoomu | Analysis of Series of Trials Where Varieties Differ in Sensitivity to Locations[END_REF]. In these models, the interaction between germplasm i and environment j was modelled as a multiplicative term η i θ j contributing to the expectation term with the remaining part adding to the residual term. The Finlay-Wilkinson coefficient is considered as both a static and a dynamic indicator of stability (Becker and Leon, 1988;[START_REF] Lin | Stability Analysis: Where Do We Stand?[END_REF]. In this model, statically stable genotypes have a coefficient close to -1. Dynamically stable genotypes have a coefficient close to zero, but having a coefficient close to zero is not sufficient to determine dynamic stability, this also depends on the amount of G × E variation that remains unexplained by the model.

Residual term

In models ADHn and FWHn, the distribution of the residual term was assumed to be normal:

ε ij ∼ N 0, σ 2 ε ,
where N 0, σ 2 ε was the normal distribution with expectation 0 and variance σ 2 ε . However, to limit the influence of extreme values on the results of the analyses, we also developed models based on Student's t distributions. Thus, in models FWHs, FWs and ADHs, the distribution of the error term was assumed to be equal to

ε ij ∼ t 0, σ 2 ε , ν ,
where t 0, σ 2 ε , ν was the Student's t distribution with dispersion parameter σ 2 ε > 0 and ν > 2 degrees of freedom. We assumed that ν > 2 to ensure that the expectation and the variance of ε ij were defined and finite. In models FWHs, FWs and ADHs, the variance of ε ij was equal to νσ 2 ε /(ν -2). The normal distribution can be considered as a t distribution with ν tending to +∞. For additive models, the residual combined the G × E contribution and errors, i.e. experimental errors and environmental heterogeneity in each trial, while for FW models, it combined the part of G × E not explained by η and errors.

Prior distribution

The statistical analysis was carried out in a Bayesian framework, so that a joint prior distribution was placed on model parameters. We placed weakly informative priors on σ ε and ν [START_REF] Cao | Bayesian Inference of Spatially Correlated Random Parameters for On-Farm Experiment[END_REF]Gelman, 2006;[START_REF] Juárez | Model-Based Clustering of Non-Gaussian Panel Data Based on Skew-t Distributions[END_REF]:

σ ε ∼ N + (0, λ 2 ε ), ν ∼ Γ(2, 0.1),
with ν > 2, and where λ ε was a known prior value of the standard deviation of the trait, N + (0, λ 2 ε ) was the normal distribution restricted to positive values with parameters 0 and λ 2 ε , and Γ(2, 0.1) was the gamma distribution with shape parameter 2 and rate parameter 0.1.

Given the high data disequilibrium and the large numbers of germplasm and environments, we decided to implement a hierarchical Bayesian approach. In all the models except the FWs model, α i , θ j and when present η i were assumed to follow hierarchical distributions:

α i ∼ N µ Y , σ 2 α , η i ∼ N 0, σ 2 η , θ j ∼ N 0, σ 2 θ ,
where µ Y , σ α , σ η and σ θ were unknown parameters. Then, we placed weakly informative prior distributions on the hyperparameters µ Y , σ α , σ η and σ θ :

µ Y ∼ N λ µ , λ 2 ε , σ α ∼ N + (0, λ 2 ε ), σ θ ∼ N + (0, λ 2 ε ), σ η ∼ N + (0, 0.75 2 ),
where λ µ was a known prior value of the trait mean. Germplasm main effects, environment main effects, germplasm sensitivities and residuals were assumed to be independent given the hyperparameters, σ ε and ν. In model FWs, the hierarchical distributions of α i , η i and θ j were replaced by weakly informative prior distributions:

α i ∼ N µ Y , λ 2 ε , η i ∼ N 0, 0.75 2 , θ j ∼ N 0, λ 2 ε .
The values chosen for λ ε and λ µ are in Appendix A.1.

Posterior distribution

Bayesian inference is based on the posterior distribution of the model parameters. This distribution was estimated using Markov chain and Monte Carlo (MCMC) methods. These methods simulate the values of the model parameters according to a Markov chain that converges to the posterior distribution of these parameters [START_REF] Robert | The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation[END_REF]. The MCMC methods were implemented using R (R Core Team, 2014) and the package rstan (Stan Developpement Team, 2016), that performs Hamiltonian Monte Carlo (HMC) sampling. This method aims at reducing the correlation between successive sampled values by using a proposal distribution based on Hamiltonian dynamics [START_REF] Neal | MCMC Using Hamiltonian Dynamics[END_REF]. Four MCMC chains were run independently to test for convergence. The initial values of each chain were taken randomly. For each chain, the burn-in consisted of 200 iterations, then 5,000 iterations were performed for all models, except FWs where 8,000 iterations were required. The average calculation time (for a given trait and a given model) was 9 minutes and the maximum time was 22 minutes (with FWs), with a computer intel CORE i7©. Estimates of the Gelman-Rubin statistic were smaller than 1.02 and the effective sample size was greater than 400 for each parameter in all tested models.

Model comparison

We compared the predictive ability of models using leave-one-out cross-validation, which seems more appropriate than Bayes factors for selecting models that approximate the process generating the data [START_REF] Lartillot | Identifying the best approximating model in Bayesian phylogenetics: Bayes factors, crossvalidation or wAIC?[END_REF]. We estimated the expected logarithmic predictive density using the R package LOO [START_REF] Vehtari | Practical Bayesian Model Evaluation Using Leave-One-out Cross-Validation and WAIC[END_REF].

This criterion was equal to

elpd loo = (i,j)∈C ln(p(Y ij |Y -ij )),
where Y -ij was the dataset without observation Y ij , and p(Y ij |Y -ij ) was the leave-one-out posterior density of Y ij . The larger this criterion, the better the agreement between the model and the data. The elpd loo criterion was also used to identify extreme observations. The quantity ln(p(Y ij |Y -ij )) can be understood as the contribution of observation Y ij to elpd loo . Observations with low contributions are unlikely and can be considered extreme observations.

For main effects and sensitivities, we estimated the average standard deviation of estimates, which allowed us to have an estimate of the precision of these effects. To be able to compare the precision between traits, for α and θ we estimated the average coefficient of variation by dividing this standard deviation by the general average µ Y .

Data analysis

Model parameters were studied using the best model as determined by the methods described above.

Variance decomposition

In order to quantify the influence of model terms on observations, the variance of an observation was decomposed. Since α i , θ j , η i and ϵ ij were assumed to be conditionally independent, the variance of an observation given the hyperparameters, σ 2 ε and ν was equal to

V ar(Y ij ) = V ar(α i + θ j + η i θ j + ε ij ) = σ 2 α + σ 2 θ + σ 2 η σ 2 θ + V ar(ε ij ).
As the best model involved the t distribution, V ar(ε ij ) was equal to νσ 2 ϵ /(ν -2). The proportions of variance explained by the germplasm main effect, the environment main effect and the interaction effect were equal to

π(α) = σ 2 α V ar(Y ij ) , π(θ) = σ 2 θ V ar(Y ij ) , π(ηθ) = σ 2 η σ 2 θ V ar(Y ij ) .
π(α) is also called broad-sense heritability. The proportion of variance explained by the model (coefficient of determination) was equal to

R 2 = π(α) + π(θ) + π(ηθ) = σ 2 α + σ 2 θ + σ 2 η σ 2 θ V ar(Y ij ) .
This definition of R 2 ensured that R 2 ≤ 1 [START_REF] Gelman | R-Squared for Bayesian Regression Models[END_REF]. We also estimated the proportion of the variance of G × E interactions and experimental errors that was explained by the η i θ j term, defined by

ρ = V ar(η i θ j ) V ar(η i θ j + ε ij ) = σ 2 η σ 2 θ σ 2 η σ 2 θ + V ar(ε ij )
.

Characterization of germplasm

Germplasm main effects and sensitivities were estimated. In addition, we estimated two stability indicators, the static stability S 2 i (Becker and Leon, 1988) and the ecovalence W i [START_REF] Wricke | Uber Eine Methode Zur Erfassung Der Okologischen Streubreite in Feldverzuchen[END_REF] which is an indicator of dynamic stability. Static stability describes the response of a genotype that maintains a constant performance across environments, while dynamic stability describes the response of a genotype showing a constant difference with the average response of all genotypes tested in each environment [START_REF] Annicchiarico | Genotype X Environment Interactions: Challenges and Opportunities for Plant Breeding and Cultivar Recommendations[END_REF]. Due to data imbalance, the empirical estimates of these indicators were biased. Thus, we defined stability indicators by means of theoretical variances from the Bayesian model described above [START_REF] Cotes | A Bayesian Approach for Assessing the Stability of Genotypes[END_REF][START_REF] Piepho | Stability Analysis Using the SAS System[END_REF].

Using the independence assumptions of the model, we obtained for germplasm i,

W i = V ar(η i θ j + ε ij ) = η 2 i σ 2 θ + V ar(ε ij ), S 2 i = V ar(θ j + η i θ j + ε ij ) = (1 + η i ) 2 σ 2 θ + V ar(ε ij ) = (1 + 2η i )σ 2 θ + W i .
The larger these indicators, the less stable the germplasm. Becker (1981) applied the same decomposition with the empirical variances. These stability indicators are approximations of the static stability and ecovalance in a balanced framework.

Moreover, we tested whether the "type" of germplasm (cross, landrace, registered variety, mixture of germplasm and historic variety) had an influence on germplasm parameters (α i , η i , S 2 i and W i ) by running a one-way ANOVA and Tukey-Kramer HSD test with germplasm type as factor.

Results

Model comparison

Predictive capacity of models

According to the elpd loo criterion, the non-hierarchical FWs model was less predictive than the hierarchical FWHs model for all the traits (Fig. 1). Using the latter model shrank the estimates of η and sometimes α (Fig. 2). With the non-hierarchical model (FWs), some estimates (α i and η i ) seemed to be unreliable, in particular some germplasm means were extreme and some FW coefficients were larger than 1 or smaller than -1.

The hierarchical models with a Student likelihood (FWHs, ADHs) were more predictive than the models with a normal distribution (FWHn, ADHn), all the more as ν was low (Tab. 5). For protein content, ν was equal to 20, so the t distribution was close to a normal distribution. The t distribution reduced the shrinkage of FW coefficients (Fig. 3). Moreover, t models better accounted for extreme data than normal models (Fig. 4). These extreme data mainly came from germplasm that were not replicated in the trials.

The Finlay-Wilkinson models (FWHs, FWHn) were slightly more predictive than the simple additive models (ADHs, ADHn), except for protein content, where the difference was not significant (Fig. 1). This difference was smaller than the differences due to the distribution of residuals and the hierarchization of parameters.

The elpd loo criterion was estimated using Pareto smoothed importance sampling [START_REF] Vehtari | Practical Bayesian Model Evaluation Using Leave-One-out Cross-Validation and WAIC[END_REF]. This method tends to be less precise for models that do not fit the data well. Thus, as expected, estimates of elpd loo were more reliable for the two hierarchical models with a t likelihood (FWHs and ADHs) than for the other models, in particular model FWs (Supplementary Tab. B.1).

Precision of estimates and distribution of residuals

For the models with a t distribution, the estimate of the number of degrees of freedom (ν) varied between 3.8 and 28.2 (close to a normal distribution) (Tab. 5). Thus, the shape of the distribution of residuals depended on the trait. This result confirmed that the number of extreme observations was not negligible in our data, and that models with a t distribution were more appropriate. In the latter case, the variation ranges of residuals were wider but with more residual values close to 0 for the t distribution than the normal distribution (Fig. 3).

Models had similar estimate precision, except for model FWs, which had less precise estimates. This result confirmed that a basic joint regression, i.e. non-hierarchical model, was not suited to our unbalanced data.

Parameters α and θ were estimated more precisely (difference in coefficient of variation between -0.1 and 0.4, Tab. 5) for t models (ADHs and FWHs) than for normal models (ADHn and FWHn). This result was consistent with Fig. 4, where extreme data were better predicted by FWHs than by FWHn, except for protein content.

-8600 -8500

-8400

Plant height

ELPD A D H n A D H s F W H n F W H s F W s -550 -500 -450 -400 Spike weight ELPD A D H n A D H s F W H n F W H s F W s -1780 -1740 -1700 -1660 Protein ELPD A D H n A D H s F W H n F W H s F W s -5500 -5400 TKW ELPD A D H n A D H s F W H n F W H s F W s Figure 1.
Predictive capacity of models. elpd loo and its associated standard error for the four studied traits.

Data analysis

In the following, we used the FWHs model which proved to be the best model in terms of prediction and accuracy of the estimated parameters.

Variance decomposition

The proportion of variance explained by each term of the model depended on the trait (Tab. 6). For all four traits, the environment effect was highly explanatory. For height and TKW, a relatively large part of the total variance was explained by the germplasm effect (resp. 23.8% and 16.1% ), whereas this part was much smaller for spike weight and protein content (10.9% and 5.6%). The proportion of variance explained by the sensitivity effect η was not significantly different from 0 for protein content and low for the three other traits. It explained 6.5%, 4.8% and 6.9% of the variance of G × E interactions and experimental errors (ρ parameter)

for plant height, spike weight and TKW, respectively.

Characterization of germplasm 336

The correlation between germplasm sensitivity (η i ) and static stability (S 2 i ) was very close to 1 for all traits 337 while germplasm sensitivity was poorly correlated to W i (Tab. 7). The main effect α i had a low correlation 338 with η i and S 2 i , except for plant height and spike weight. Correlations between W i and α i were low and in 339 most cases not significant.

340

Plant height was found to depend on the type of germplasm, landraces being taller than historic varieties, 341 which were themselves taller than registered varieties. Registered varieties were significantly more stable 342

(static stability and FW coefficient) than everything other than mixtures, but less stable dynamically (ecov-343 alance). No germplasm parameters were significantly dependent on germplasm type for protein content and 344 spike weight. TKW germplasm main effects did not depend on the type of germplasm, but for this trait lan-345 draces and mixtures appeared statically more stable than historic varieties.

Discussion

347

To fit the characteristics of PPB trials, i.e., few inter-farm replicates and possible extreme data, we de-348 veloped several models and we found that the hierarchical Finlay-Wilkinson model with t residuals was the 349 best for prediction and parameter precision. Then we compared the performance and stability of different 350 germplasm types. to the elpd loo criterion. Black (resp. red) dots correspond to observations that were measured on germplasm that were replicated (resp. not replicated) within trials.

addition, its estimates had poor precision and it led to extreme sensitivity estimates, with values close to 1 or -1 (Fig. 2).

In contrast, hierarchical joint regression appeared more suited to our data structure. Model FWHs had the largest elpd loo values for three traits out of four. Placing a hierarchical distribution on sensitivities constrained estimates and brought them closer to 0. This led to more satisfactory sensitivity estimates, since they were well below 1 in absolute value.

Three strategies have previously been used to manage incomplete G × E data: i) subset the total dataset to obtain an almost balanced subset for the analysis [START_REF] Ceccarelli | Decentralized-Participatory Plant Breeding: An Example of Demand Driven Research[END_REF], ii) predict missing data with a more or less complex model and use these predictions in the analysis [START_REF] Kumar | High-Yielding, Drought-Tolerant, Stable Rice Genotypes for the Shallow Rainfed Lowland Drought-Prone Ecosystem[END_REF][START_REF] Woyann | Estimation of Missing Values Affects Important Aspects of GGE Biplot Analysis[END_REF], and iii) use a model more robust to unbalanced data, provided it complies with model validation conditions [START_REF] Assis Tog De | A Weighted AMMI Algorithm for Nonreplicated Data[END_REF][START_REF] Van Frank | Influence of Experimental Design on Decentralized, on-Farm Evaluation of Populations: A Simulation Study[END_REF]. We used the last strategy to maximise the amount of information from the data (less data excluded than in the first strategy) with a one-step process (unlike the second strategy). Estimates were made using model FWHs and are given in %. Mean: mean of the posterior distribution; 95% CI: 95% credible intervals. R 2 is the coefficient of determination. π(α), π(θ) and π(ηθ) are respectively the proportion of variance explained by α, θ and ηθ. ρ is the proportion of the variance of G × E and errors explained by ηθ.

van Frank et al. ( 2019) investigated the influence of MET design on the evaluation of germplasm using model FWHn and simulations. They found that when data were highly unbalanced, this evaluation was more reliable with at least 100 environments and when some germplasm were replicated over at least 5 environments per germplasm. Therefore, we restricted the analysis to a subset of data that met these conditions. Having 100 environments seems difficult to achieve, but it should be noted that this number combines the number of trials per year and the number of years, which may be achieved by a network of many farmers experimenting over several years. [START_REF] Cotes | A Bayesian Approach for Assessing the Stability of Genotypes[END_REF] used a Bayesian approach to estimate FW coefficients in a MET study in order to take prior information on germplasm coming from other studies into account. A similar approach was used by [START_REF] Couto | Eberhart and Russel's Bayesian Method in the Selection of Popcorn Cultivars[END_REF], [START_REF] Foucteau | Statistical Analysis of Successive Series of Experiments in Plant Breeding: A Bayesian Approach. Quantitative genetics and breeding methods: the way ahead[END_REF][START_REF] Foucteau | Statistical Analysis of Successive Series of Experiments in Plant Breeding: A Bayesian Approach. Quantitative genetics and breeding methods: the way ahead[END_REF][START_REF] Nascimento | Bayesian Segmented Regression Model for Adaptability and Stability Evaluation of Cotton Genotypes[END_REF] and was found to greatly improve the results. Here, we used little prior information. But in the future, previous evaluation studies may provide n: number of germplasm, α i : mean germplasm effect, η i :mean sensitivity (FW coefficient), S 2 i : mean static stability and W i : mean ecovalence.

stronger prior information on germplasm behaviour.

Extreme observations

Extreme observations were more frequent in our dataset than expected under the normal distribution for three traits out of four (Fig. 4). For these traits, using a t distribution increased elpd loo values, and the estimate of the number of degrees of freedom of this distribution was smaller than 10 (Tab. 5). These extreme observations could occur in our dataset for several reasons: because most of the populations were not replicated within the trials, because cultivation environments were less controlled, or because a non-negligible part of the GxE interaction was not captured by the multiplicative term of the FW model. The normal distribution was appropriate for the trait protein content. It is difficult to explain why this trait had fewer extreme observations. A possible explanation could be that the measurement of protein content is more standardized than other trait measurements. For plant height, extreme values occurred only for non-replicated micro-plots with a global measurement and never with data from the average of 25 plants (Sect. 2.2.2), suggesting that this measurement is less accurate. For TKW, the kernel count could be affected by broken kernels due to over-drying or incorrect threshing settings leading to an overestimation of the number of kernels in the sample. Another possible explanation is that protein content is less variable under different conditions than plant height and spike weight [START_REF] Kazakou | Are Trait-based Species Rankings Consistent across Data Sets and Spatial Scales[END_REF].

Using a t distribution did not affect the estimates of germplasm and environment main effects. On the contrary, it improved the estimates of sensitivities. It reduced their shrinkage and allowed the multiplicative term of the FW model to better capture G × E interactions (Fig. 3).

The Student distribution is expected to take better account of extreme data and to yield more robust estimates (Besag and Higdon, 1999;[START_REF] Lange | Robust Statistical Modeling Using the t Distribution[END_REF][START_REF] Rosa | Robust Linear Mixed Models with Normal/Independent Distributions and Bayesian MCMC Implementation[END_REF]. Extreme data are more likely to occur when varieties are not replicated within trials, which is frequent in this dataset (Fig. 4). [START_REF] Rosa | Robust Linear Mixed Models with Normal/Independent Distributions and Bayesian MCMC Implementation[END_REF] found that a normal likelihood underestimated a sex effect compared to a t likelihood. This effect was estimated less precisely with a normal distribution, which is consistent with our results for plant height, spike weight and TKW. A Student likelihood appears to be a good solution for dealing with extreme data, in particular in stability analyses, where extreme observations are sometimes removed (this is justified when they are extreme because of experimental errors, but not when they are due to natural variability). While this distribution has recently been used to implement robust alternatives to BLUP [START_REF] Gianola | Prediction of complex traits: robust alternatives to best linear unbiased prediction[END_REF] or to handle environmental heterogeneity in a single trial [START_REF] Cao | Bayesian Inference of Spatially Correlated Random Parameters for On-Farm Experiment[END_REF], to our knowledge, it has not already been used in MET studies.

Computing time

Series of trials often include many genotypes and environments, leading to large data sets. Thus, their analysis using mixed or hierarchical models is generally computationally demanding [START_REF] Smith | The Analysis of Crop Cultivar Breeding and Evaluation Trials: An Overview of Current Mixed Model Approaches[END_REF]. The computational load can be reduced by using approximate estimation methods [START_REF] Nabugoomu | Analysis of Series of Trials Where Varieties Differ in Sensitivity to Locations[END_REF] or efficient algorithms, such as algorithms based on sparse matrix operations [START_REF] Gilmour | Average Information REML: An Efficient Algorithm for Variance Parameter Estimation in Linear Mixed Models[END_REF][START_REF] Thompson | A Sparse Implementation of the Average Information Algorithm for Factor Analytic and Reduced Rank Variance Models[END_REF]. Hierarchical joint regression has already been implemented using Gibbs sampling or Jags [START_REF] Lian | FW: An R Package for Finlay-Wilkinson Regression That Incorporates Genomic/Pedigree Information and Covariance Structures Between Environments[END_REF][START_REF] Van Frank | Genetic Diversity and Stability of Performance of Wheat Population Varieties Developed by Participatory Breeding[END_REF]. Our implementation based on Hamiltonian Monte Carlo and Stan was more efficient since it required fewer iterations [START_REF] Gelman | Analysis of Variance-Why It Is More Important than Ever[END_REF]. It allowed us to analyze large datasets in about 10 minutes.

To reduce computing time, the analyses were carried out in two steps. First, germplasm means were estimated using within-trial analyses. Then, these estimates were gathered and analyzed using a between-trial analysis. Thus, this two-stage approach analyzes G × E means without taking account of their standard error. Rivière et al. (2015a) developed a flexible method for estimating the experimental variance of trials with low intra-farm replication (farm design presented in Tab. 2). An easy way to integrate the variability estimated in the first stage would be for instance to use the same method as [START_REF] Couto | Eberhart and Russel's Bayesian Method in the Selection of Popcorn Cultivars[END_REF].

Main effects

For the four traits studied, we found that the environmental part of the variance was large (from 51% for TKW to 77% for protein content, Tab. 6), which is consistent with the diversity of the cropping environments encountered (soil, climate, cropping practices...). Nevertheless, heritability was still significant with plant height > TKW > spike weight > protein. Rivière et al. (2015b) found (with data included in our study) a similar ranking in heritability: plant height > TKW = protein > spike weight. Plant height is known to be quite heritable due to a relatively simple genetic architecture with a few major genes, such as the well known Green Revolution Rht1 and Rht2 genes [START_REF] Peng | Green Revolution'Genes Encode Mutant Gibberellin Response Modulators[END_REF]. In our study, the presence of both quite recent registered varieties and varieties dating from before the second World War, very likely led to varieties containing different alleles for these loci and increased variability for height. The decrease in plant height from landraces to historic varieties and registered varieties appears very clearly (Tab. 8) as also found in several studies (Bektas et al., 2016;[START_REF] Cantarel | Plant Functional Trait Variability and Trait Syndromes among Wheat Varieties: The Footprint of Artificial Selection[END_REF].

Germplasm stabilities

FW coefficients explained a low proportion of the total variance (between 0.2% and 2.2%) and a low proportion of the variance of G × E interactions and errors (between 1.2% and 6.9%, Tab. 6). We can presume that the explanation of the interaction by the FW parameter is weaker the greater the number of environments for example, 29% for 12 studies with less than 10 environments, and 12% for studies with more than 10 environments (Brancourt-Hulmel et al., 1997). Other classical models, such as AMMI (additive main effect and multiplicative interaction) or GGE (G + G × E) models, might explain a larger part of GxE interactions.

Missing data estimation methods allow these models to be used when the data is highly unbalanced, with up to from 40% unbalanced data for a MET with less than 20 environments to 60% unbalanced data for MET with at least 40 environments [START_REF] Woyann | Estimation of Missing Values Affects Important Aspects of GGE Biplot Analysis[END_REF][START_REF] Yan | Biplot Analysis of Incomplete Two-way Data[END_REF]. However, these datasets are more balanced than ours, and, as found by [START_REF] Rodrigues | A Comparison between Joint Regression Analysis and the Additive Main and Multiplicative Interaction Model: The Robustness with Increasing Amounts of Missing Data[END_REF], FW is more robust than AMMI when the data are highly unbalanced (75%). In our study, most germplasm occurred in a limited number of environments, so that a parsimonious and very simple modelling of G×E interactions had to be used. An alternative approach would be to better characterize the environments and thus explain the environmental effects and part of the G × E interaction using environmental variables (Piepho, 2022).

Although sensitivities explained a rather low proportion of variance, HFWs model had larger elpd loo values than additive models for three traits out of four. In addition, for these traits, some sensitivity estimates were not negligible, with values close to 0.2 or 0.3. Interaction effects then represented 20% or 30% of environmental effects. Additive models were appropriate for the protein content trait. It was found that the multiplicative term of the FW model was not significant for protein content, both in a balanced network of 15 environments in Serbia [START_REF] Hristov | Genotype by Environment Interactions in Wheat Quality Breeding Programs in Southeast Europe[END_REF] and in 12 environments in Swiss organic trials [START_REF] Knapp | Quality Traits in Winter Wheat: Comparison of Stability Parameters and Correlations between Traits Regarding Their Stability[END_REF]. On the contrary, [START_REF] Mut | Stability of Some Quality Traits in Bread Wheat (Triticum Aestivum) Genotypes[END_REF] found significant FW coefficients for a balanced network of 7 environments in Turkey.

These contrasting results could be explained by differences between numbers of environments or between genetic diversities.

For plant height, we found that registered varieties were more statically stable but less dynamically stable (Tab.8). This can be explained by the fact that there are only a few registered varieties in the trials, therefore they have little influence on the average height, which can fluctuate greatly between trials, and therefore the deviation from this average will be greater for this type.

Static and dynamic stabilities were difficult to estimate since our series of trials was very unbalanced. In particular, raw estimates of these stabilities were not reliable, since they were much influenced by the unbalanced nature of the data. By using theoretical variances, the FW model allowed us to calculate simple indicators of static and dynamic stability in the wheat PPB dataset. To our knowledge, the FW model has never been used for this purpose before.

Dependence between stability and mean is widespread [START_REF] Reckling | Methods of Yield Stability Analysis in Long-Term Field Experiments[END_REF], but in our case, the correlation was low, which simplified interpretation of the stability analysis. Several studies for different traits and with balanced MET found a very strong correlation between FW coefficient and the static stability (Becker, 1981;[START_REF] Fasahat | An Overview on the Use of Stability Parameters in Plant Breeding[END_REF][START_REF] Reckling | Methods of Yield Stability Analysis in Long-Term Field Experiments[END_REF]. However, in our case, this relationship was even stronger (Tab. 7), probably because of the assumption that the variance of residuals did not depend on the genotype.

As in many other studies, the residual variance was assumed to be independent of germplasm throughout our study. Allowing the residual variance to depend on the genotype could improve the estimates of stability indicators [START_REF] Cotes | A Bayesian Approach for Assessing the Stability of Genotypes[END_REF][START_REF] Couto | Eberhart and Russel's Bayesian Method in the Selection of Popcorn Cultivars[END_REF]. In particular, the dynamic indicator would be similar to the Shukla Stability Variance, i.e, the varietal variance of G × E interactions [START_REF] Cotes | A Bayesian Approach for Assessing the Stability of Genotypes[END_REF]. However, estimating a residual variance and a FW coefficient for each germplasm could be difficult in our study, as most of the germplasm appeared in only a few environments.

When relationships were significant, mixtures were always in a more stable (statically and dynamically) statistical group (Tab. 7). This result supports the fact that within-plot diversity stabilizes performances [START_REF] Döring | Comparative Analysis of Performance and Stability among Composite Cross Populations, Variety Mixtures and Pure Lines of Winter Wheat in Organic and Conventional Cropping Systems[END_REF][START_REF] Kiaer | Effects of Inter-Varietal Diversity, Biotic Stresses and Environmental Productivity on Grain Yield of Spring Barley Variety Mixtures[END_REF].

In the wheat PPB program, the populations tested were heterogeneous and their genotypic composition could vary over years and farms [START_REF] David | Bayesian Inference of Natural Selection from Spatiotemporal Phenotypic Data[END_REF]. In this analysis, such variations were considered as part of the response of a population to a given environment for the sake of simplicity. Therefore the G × E interactions could be overestimated (resp. underestimated) if populations underwent diversifying (resp.

stabilizing) selection pressures within farms.

Role of statistical methods in the wheat PPB project

This work was developed following the co-construction of an experimental set-up suitable for decentralized on-farm evaluation and selection, and research into the best methods for analyzing the resulting data. It was part of the methodology we set up in a wheat PPB program, which was based on a collaboration between farmers, associations of farmers and researchers [START_REF] Dawson | Collaborative Plant Breeding for Organic Agricultural Systems in Developed Countries[END_REF]. The farmers could freely choose the populations they wanted to test, so that a wide genetic diversity could be evaluated in a wide range of environments. There were on average more than 130 environments resulting from the combination of years and farms. The number of genotypes evaluated was large compared to other studies, but it was smaller than in CIMMYT's MET, which involved between 500 and 800 genotypes tested in 12 MET between 1945[START_REF] Lin | Stability Analysis: Where Do We Stand?[END_REF](Braun et al., 1997).

One aim of the project was to provide farmers with information to help them select new germplasm for testing on their farm. The statistical tools we developed sought to cope with the large degree to which this series of trials was unbalanced. Their objectives were the same as in other MET analyses : (i) estimate and predict germplasm' values for traits of interest for breeding, (ii) study the stability of germplasm over several environments, (iii) select new germplasm to be tested in new locations [START_REF] Cotes | A Bayesian Approach for Assessing the Stability of Genotypes[END_REF]. MET are usually carried out to find stable germplasm that perform well on average over many locations, or to detect special local adaptations to certain environments [START_REF] Annicchiarico | Defining Subregions and Estimating Benefits for a Specific-Adaptation Strategy by Breeding Programs: A Case Study[END_REF][START_REF] Gauch | Statistical Analysis of Yield Trials by AMMI and GGE: Further Considerations[END_REF]. Here, while farmers were mostly interested in selecting the best germplasm adapted to their local pedo-climatic conditions, farming practices and marketing objectives, information retrieved from the farmers' network on new varieties to introduce in their trials could also be useful.

Conclusion

The proposed hierarchical model was aimed at improving the estimates of the parameters of the FW model from unbalanced datasets. This model was complex and was easier to implement in a Bayesian framework.

Placing hierarchical distributions on model parameters and modelling residuals using a t distribution improved the estimates of main and interaction effects. This model allowed us to estimate static and dynamic stability indicators despite the high level of unbalanced data. Main effects and stability indicators provide information on the behaviour of genotypes in different environments, which farmers could use in their selection process.

Participatory research raises new research questions and contributes to the development of new methods for societal action [START_REF] Kastenhofer | Sustaining Sustainability Science: The Role of Established Inter-Disciplines[END_REF]. In PPB programs, all the methodology is based on collective and collaborative work and action between farmers, associations of farmers and researchers (Brac de la Perrière et al., 2011). New statistical methods can contribute to a better use of such complex multi-environment data in the selection process, and more generally to the effectiveness of participatory research (Martin and [START_REF] Sherington | Participatory Research Methods-Implementation, Effectiveness and Institutional Context[END_REF].

Figure 2 .

 2 Figure 2. The first column presents the distribution of the trait to be explained (in grey). The last four columns compare the hierarchical (red) and non hierarchical (blue) versions of the FW model with a Student law for the residuals, and show the smoothed histograms of main effects and FW coefficients.

Figure 3 .

 3 Figure 3. Comparison of hierarchical FW models with different residual laws, the Student (red) and the normal (blue). These graphics show the smoothed histograms of main effects and FW coefficients.

-

  the data from a highly unbalanced series of trials352As the farmers of the program chose the germplasm they assessed, the data obtained from the series of 353 trials were very unbalanced, with more than 90% of the G×E combinations missing. This made the estimation 354 of germplasm main effects and sensitivities difficult. Although the Finlay-Wilkinson model was parsimonious, 355 a basic joint regression with weakly-informative prior distributions (model FWs) was not able to cope with 356 this level of disequilibrium. According to the elpd loo criterion, model FWs was not the best model (Fig.1). In 357

Figure 4 .

 4 Figure 4. Comparison of t and normal models (FWHs vs FWHn) in terms of the contributions of observations

  Figure B.1. Comparison of models FWHs and ADHs for the distribution of germplasm main effects (α), environment main effects (θ), FW coefficients (η) and residuals (ε) for each trait. Red: model FWHs; Blue: model ADHs.

Figure B. 2 .

 2 Figure B.2. Comparison of models ADHs and ADHn in terms of the contributions of observations to the elpd loo criterion. Black (resp. red) dots correspond to observations that were measured on germplasm that were repeated (resp. not repeated) within trials

Table 1 .

 1 Main notations.

	Notation Meaning
	PPB	Participatory plant breeding
	OA	Organic agriculture
	RSP	French farmers' seed network
	MET	Multi-environment trial
	G × E	Genotype × environment interaction
	FW	Finlay Wilkinson
	MCMC	Markov chain Monte Carlo
	LOO	Leave one out
	α	Germplasm main effect
	θ	Environment main effect
	η	Germplasm sensitivity (FW coefficient)
	S 2	Germplasm static stability
	W	Germplasm ecovalence (a dynamic stability)
	elpd loo	LOO expected logarithmic predictive density

Table 2 .

 2 Experimental designs of the 189 trials used in the statistical analysis. Nb: number, Environment:

			Nb of repeated	Nb of gemplasms by	Nb of environ-
			germplams	environment	ments
	Complete blocks	2 to 3	6 to 45	7 to 45	24
	Incomplete blocks	3 to 4	3 to 49	6 to 58	11
	Regional farm	2	5 to 16	7 to 81	20
	Statellite farm	1	1 to 22	5 to 79	102
	Unreplicated	1	0	5 to 39	32

Table 3 .

 3 Description of the dataset. Disequilibrium: proportion of missing values in the Germplasm x Environment table in %.

	Trait	Observations Germplasm Environments Disequilibrium Farms Years
	Plant height	1437	124	117	90	44	11
	Spike weight	1804	172	148	93	52	10
	Protein	1332	144	111	92	44	9
	TKW	1982	177	165	93	58	11

Table 4 .

 4 The five models fitted.

		Expectation term Residual term Prior distribution
	ADHn	Additive	Normal	Hierarchical
	ADHs	Additive	Student	Hierarchical
	FWHn Finlay Wilkinson	Normal	Hierarchical
	FWHs	Finlay Wilkinson	Student	Hierarchical
	FWs	Finlay Wilkinson	Student	Weakly informative

Table 5 .

 5 Number of degrees of freedom and precision of estimates.

		Trait	Model	ν	cv(α) cv(θ) sd(η)
			ADHn		3.0	2.8
			ADHs	3.8 (0.5)	2.8	2.6
		Plant height	FWHn		3.1	2.9	0.08
			FWHs	3.5 (0.4)	2.8	2.7	0.09
			FWs	3.3 (0.4)	4.7	4.5	0.20
			ADHn		5.3	5.2
			ADHs	8.1 (2.2)	5.2	5.0
		Spike weight FWHn		5.3	5.2	0.12
			FWHs	8 (2.2)	5.3	5.1	0.11
			FWs	10.2 (4.1)	6.8	5.7	0.29
			ADHn		2.7	2.7
			ADHs	19.9 (9.6)	2.6	2.7
		Protein	FWHn		2.6	2.7	0.05
			FWHs	19.6 (9.6)	2.7	2.7	0.05
			FWs	28.2 (13.4)	4.8	4.1	0.25
			ADHn		2.8	2.8
			ADHs	4.2 (0.5)	2.6	2.5
		TKW	FWHn		2.8	2.8	0.15
			FWHs	4 (0.5)	2.7	2.5	0.17
			FWs	3.8 (0.4)	3.5	3.2	0.33
		Plant height	Spike weight		Protein	TKW
		Mean 95% CI	Mean 95% CI	Mean 95% CI	Mean 95% CI
	R 2	87.1 82.9∼90.2	78.2 73.8∼82.1	83 79.1∼86.7	69.7 64.1∼74.8
	π(α)	23.8 17.9∼30.7	10.9 7.7∼14.7	5.6 3.7∼8.2	16.1 12∼20.8
	π(θ)	62.4 54.5∼69.8	66.2 60.1∼71.9	77.2 72∼82	51.4 44.6∼58.1
	π(ηθ)	0.9 0.4∼1.6	1.1 0.4∼2.1	0.2 0∼0.8	2.2 1∼3.8
	ρ	6.5 2.8∼12	4.8 1.7∼9.3	1.2 0∼4.7	6.9 3.1∼12

ν: posterior means, with posterior standard deviations in parentheses, of the number of degrees of freedom of the t distribution; cv(α), cv(θ): average posterior coefficients of variation of germplasm and environment main effects; sd(η): average posterior standard deviation of germplasm sensitivities (FW coefficients).

Table 6 .

 6 Variance decomposition.

Table 7 .

 7 Correlation between germplasm parameters. *, **, *** : significant at P = 0.05, P = 0.01, P = 0.001 respectively. α i : germplasm effect, η i : germplasm sensitivity (FW coefficient), S 2 i : static stability, W i : ecovalence.

			Trait	Pearson correlation between
				αi|ηi	αi|S 2 i	αi|Wi ηi|S 2 i	ηi|Wi S 2 i |Wi
		!	Plant height Spike weight 0.34*** 0.35*** 0.24* 0.43*** 0.4*** -0.41	0.997*** -0.34 0.999*** 0.17	-0.26 0.21
			Protein	0.13	0.14*	-0.01	0.999*** -0.01	0
			TKW	0.23*	0.23*	0.12	0.996*** 0.27	0.36
	Trait		Registered		Historic	Landrace	Cross	Mixture	P-value
		αi	862 c		1136 b	1221 a	1175 ab	1189 ab	<0.001
	Plant height	ηi	n=6 -0.11 b	n=16 -0.01 a	n=21 0 a		n=74 0.01 a	n=7 -0.01 ab	0.001
		S 2 i	37566 b		44170 a	44616 a	45556 a	43911 ab	0.005
		Wi	8835 a		7680 b	7678 b	7716 b	7707 b	<0.001
		αi	2.03		1.99	1.94	1.96	2.02	0.439
	Spike weight ηi	n=8 0.02	n=20 0.02	n=40 -0.01	n=89 0	n=15 -0.01	0.297
		S 2 i	0.34		0.34	0.32	0.33	0.32	0.303
		Wi	0.08		0.08	0.08	0.08	0.08	0.712
		αi	11.1		11.4	11.4	11.4	11.5	0.125
	Protein	ηi	n=9 0.001	n=27 0.001	n=34 -0.001	n=62 0	n=12 0.001	0.935
		S 2 i	3.41		3.41	3.39	3.4	3.41	0.951
		Wi	0.61		0.61	0.61	0.61	0.61	0.641
		αi	43.6		43.9	43.1	43.4	43.9	0.681
	TKW	ηi	n=8 -0.03 ab	n=20 0.08 a	n=42 -0.03 b	n=92 0.01 ab	n=15 -0.05 b	0.003
		S 2 i	32.9 ab		37.4 a	32.7 b	34.3 ab	31.7 b	0.002
		Wi	12.7		13	12.8	12.7	12.6	0.081

Table 8 .

 8 Performance and stability of types of germplasm.

P-values in bold indicate significant differences (P<0.05) between germplasm types for a given trait and different letters indicate significant differences (P<0.05) based on Tukey-Kramer HSD tests.

Acknowledgements

We thank Estelle Serpolet, Nathalie Galic and Sophie Pin for their great help in the measurements on participating farms and research stations. We thank all the farmers and facilitators participating in the project.

We thank Gérard Branlard from INRAE Clermont-Ferrand for his time when carrying out NIRS analysis and Jean-Marc Le Goff for this feedback.

Fundings

M. Turbet Delof was funded by Program PPR-CPA MoBiDiv (2021-2026) under grant agreement ANR-20-PCPA-0006 and INRAE (program on organic scaling METABIO and Biologie et Amélioration des Plantes departement).

Conflict of interest disclosure

The authors of this article declare that they have no financial conflict of interest with the content of this article. XXX and XXX are recommenders for PCI XXX.

Data, script and code availability

Data and script for models are available online: DOI 10.57745/SUTZ9U https://doi.org/10.57745/SUTZ9U .

Supplementary information

A Models

B Model comparison

Tab. B.1 provides supplementary information on the estimation of the elpd loo criterion.

Trait

Model k < 0.5 0.5 < k < 0.7 0.7 < k < 1 k > 1


Tobias Scheer

9 CəCj in French

19.1 Charette's analysis of CəCj [START_REF] Charette | Empty and pseudo-empty categories[END_REF]Charette ( :171-173, 2003:471-475) :471-475) wonders why schwa may be left unpronounced in [CəCj] cimetière (cim'tjère) 'cemetery', but not in [CəCj] atelier (*at'ljer) 'workshop' (this pattern is also analysed by Cavirani in Chapter 15). 1 On pp. 173/473f of these articles, she ascribes the contrast to the type of CC cluster created if schwa were left out: schwa may only drop if this cluster is a good coda-onset sequence in Standard Government Phonology [START_REF] Kaye | Constituent structure and government in phonology[END_REF][START_REF] Charette | Conditions on Phonological Government[END_REF]. This amounts to any cluster which exists in intervocalic position and is different from a branching onset (obstruent-liquid TR) and a bogus cluster (tl, dl). In the [CəCj] pattern, though, the cluster created by the non-pronunciation of schwa ([mt] in cim'tière) cannot be a coda-onset sequence since it encloses a schwa in the lexicon and resyllabification is prohibited in Government Phonology (Kaye et al. 1990:221). Therefore, Charette argues, the CCs at hand are in fact onset-onset clusters that afford a right-headed onset-to-onset government whose definition is the same as the one that holds within coda-onset sequences: any intervocalic cluster that exists in the language qualifies except TR and bogus [tl, dl].

(1) CəCj [START_REF] Charette | Empty and pseudo-empty categories[END_REF][START_REF] Charette | Empty and pseudo-empty categories[END_REF] a. cimetière b. atelier

IO-Ggvt PG PG O N O N O N O N O N O N O N O N O N | | | | | | | | | | | | | c i m ə t I ɛ r a t ə l I e
Thus in cim'tière under (1a), the [m] and the [t] are onsets and [t] governs [m], thereby circumscribing the empty nucleus enclosed which for that reason is rightfully empty: [møtøjV] is well-formed since both empty nuclei are taken care of, the first by onset-onset government, the second by government from the following V. In *at'lier under (1b) on the other hand, being bogus, [t] and [l] are unable to establish interonset government and therefore schwa cannot drop when followed by another empty nucleus: *[tøløjV] is ill-formed because the first empty nucleus is neither governed nor enclosed in a governing domain (while the second empty nucleus is governed by the following V as before).

Charette's analysis of [CəCj] is thus based on the management of the CC cluster that is left when schwa remains unpronounced, knowing that it encloses an empty nucleus and is followed by another empty nucleus in all cases: in Charette's analysis, [Cj] identifies as /[Cøj/] whose empty nucleus is governed. 2

Charette's earlier analysis of CəCj

Both data and analysis presented in the previous section have changed since Charette (1991:115f) has first studied [CəCj]. In 1991, she compares Quebec and Parisian French, reporting that while schwa is always realised in the former (cim[ə]tière, at[ə]lier), it may remain (or actually is) unpronounced in the latter (cim'tière, at'lier). That is, Charette does not make any difference between the cimetière type and the atelier type.

The contrasting behaviour of the two types had been observed by 1991, though, but as far as I can see only by Dell (1973:262f) for Parisian French: schwa may be left out in the cimetière type, but not in the atelier type. In her 1998 and 2003 articles, Charette acknowledges this difference also for Quebec French: she now reports that schwa may or may not be realised in the cimetière type, while (as in Parisian French) it is mandatory in the atelier type.

Regarding the analysis, in 1991 Charette argues that schwa cannot drop in Quebec French because the following Cj is a branching onset 2a, but that it may be left out in Paris French since the yod belongs to a complex (light) nucleus 2b (see Pöchtrager's Chapter 12 on light diphthongs in general and in French).

(2) CəCj (Charette 1991:115f) a. Quebec French (schwa mandatory) b. Parisian French (schwa may be left unpronounced)

PG PG O N O N O N O N O N O N O N O N | | | | | | | | | | | | | | x x x x x x x x x x x x x x x x x x | | | | | | | | | | | | | | | c i m ə t I ɛ r c i m t I ɛ r a t ə l I e a t l I e
Under (2a), the nucleus containing schwa cannot be reached by government from the following nucleus since in Standard GP government cannot jump over a branching onset (Charette 1991:113ff). Therefore, the ungoverned nucleus must be realised and schwa is pronounced. By contrast under (2b), only a single consonant separates the schwa-containing and the following nucleus: the yod is part of the latter. Thus nothing prevents government to fall on the nucleus containing schwa, which therefore may be left unpronounced.

As was mentioned, this analysis is based on incorrect data: in both varieties, there is a crucial difference between the cimetière type and the atelier type of words. But the analytic options are interesting: in [CəCj], schwa is mandatory when the yod belongs to a branching onset, while it may be left out in case the yod belongs to a complex nucleus. This contrasts with the analysis discussed in the previous section where [Cj] identifies as two independent onsets. Also note that in Charette's 1991 analysis, the syllabic status of the [Cj] cluster is responsible for the behaviour of schwa, while in her 1998/2003 analysis the syllabic identity of the CC cluster drives schwa deletion. It will be shown in Section 19.4 that the former is correct.

Syllabic identity of Cj

Complex nucleus

The complex nucleus analysis is a consequence of the incorrect data that Charette's 1991 analysis was based on: it accounts for the pattern where schwa in [CəCj] can be dropped no matter what (which is what Charette thought occurs in Paris). It may therefore be discarded.

This is particularly clear when considering cases such as cafet-ière (schwa may be dropped) and hôtel-ière (schwa must be realised) where the yod and the following vowel belong to a suffix (the morphological structure of cimetière and atelier is less clear): the syllabic identity of the suffix-initial [jV] sequence must be the same in both types of words since the

[jV] belongs to the same morpheme. If it identified as a complex nucleus, it should be able to govern the preceding schwa in both types of wordsbut this is not the case.

Branching onset

The branching onset analysis of [Cj] is impossible in Standard GP where resyllabification is prohibited (Kaye et al. 1990:221;Charette 1991:222, note 1): whatever the lexical identity of the yod in the suffix -ière, it would need to become the second element of a branching onset, whose first element is the root-final consonant, a simple onset in the lexicon.

In other words, a branching onset made of items from two different morphemes is impossible in Standard GP.

In this context, let us consider the parallel between [əTRV] and [əCjV] that Charette (1991:116f) makes explicit. In her 1991 analysis, the reason why in Quebec French schwa cannot drop in la s[ə]crétaire 'the secretary' (əTRV) and at[ə]lier (əCjV) is the same: its nucleus cannot be reached by government, which is blocked by the following branching onset (TR in the former, Cj in the latter case) (see 2a).

In French, word-internal clusters of three consonants always involve a coda followed by a branching onset [C.TR] (perdre 'to lose', arbre 'tree', filtre 'filter', mercredi 'Wednesday ', etc., Dell 1995:10-17). Therefore, knowing that the yod is indeed a consonant (rather than part of a complex nucleus, see the previous section), the [Cj] [ml+j] (sommelier). Given her analysis, Charette must say that [ʃl] and [ml] are not possible coda-onset (and thus interonset) sequences. She does not discuss the issue, but supporting evidence is the fact that these clusters do not occur in intervocalic position in French. 4 Looking at a broader record of words instantiating [CəCj], it appears that the CC cluster created by the putative loss of schwa does not bear any responsibility for the (non-)pronunciation of schwa. Rather, the pattern is governed by the consonant following schwa: if this consonant is the lateral [l], schwa cannot be left out. If this consonant is different from the lateral, schwa may be dropped.

This is shown under (3), where [CəCj] clusters produced by the suffix -ier/-ière appear. 5 The lexical record for [CəC]-ier/-ière should be near exhaustive and may contain a number of words that natives are unfamiliar with. In case speakers encounter a word that they do not know, they will treat it as a new lexical item, i.e. as a nonce word, and that does not appear to alter the grammaticality judgements reported below.

( were not included, the latter because of possible [s+C] effects). Participants had one vote for each word where schwa was left out, to be chosen among "possible", "bizarre but possible" and "impossible". With seven participants, for any given word each of the three possible votes scores between zero and seven. The average score of "possible" for [CC+j] was 4.88, against 0.4 for [Cl+j]. Conversely, "impossible" scored 0.94 with [CC+j], against 5.27 with [Cl+j] (scores for "bizarre" were comparable: [CC+j] at 1.19, [Cl+j] at 1.33). The contrast between

[CC+j] and [Cl+j] is thus very strong, in the direction predicted.

The same contrast appears when looking at percentages: 69.6% of responses for [CC+j] items were "possible" (156 of 224), against only 5.7% for [Cl+j] items (6 of 105). This difference is of course highly significant (χ 2 (2, N=224)=116,89, p<.00001)). Conversely, 13.4% of votes were on "impossible" for [CC+j], against 75.2% for [Cl+j]. The difference is also highly significant (χ 2 (2, N=224)=123,42, p<.00001).

The blame is on the lateral

The testimony of words suffixed by -ier/-ière under (3) shows that it is not the cluster type preceding yod that governs schwa loss in [CəCj] (as held by [START_REF] Charette | Empty and pseudo-empty categories[END_REF][START_REF] Charette | Empty and pseudo-empty categories[END_REF], but rather the [Cj] cluster following schwa (as argued by Charette 1991): schwa may be dropped when followed by any consonant plus yod, except if this consonant is the lateral.

Note that [START_REF] Charette | Empty and pseudo-empty categories[END_REF]Charette's ( :173, 2003:473f) :473f) generalisation according to which schwa can be dropped if the CC in [CəCj] makes a good coda-onset sequence is also contradicted by [rl] (bourrelier) and [sl] (vaisselier) where schwa is mandatory, although the clusters at hand are good coda-onset sequences in French, as shown by merle 'blackbird', arlequin 'harlequin', (avoir la) berlue 'being blind ', dirlo 'director', horloge 'clock', hurler 'to yell' etc. for [rl] and islam 'Islam', Islande 'Iceland', législatif 'legislative' etc. for [sl].

Given the parallel with [CəTR] where TR is a good branching onset (see Section 19.3.2), the generalisation suggests that yod can engage with any preceding consonant to form a branching onset, except with the lateral: speakers who can drop schwa in [CəTR] (la s'crétaire 'the secretary') can also leave it out in [C1əC2j], provided that C2 is not the lateral (cim'tière).

In other words, any [Cj] sequence is a good branching onset, except [lj] which is not. [START_REF] Charette | Conditions on Phonological Government[END_REF] was thus right in identifying the syllabic status of [Cj] in [CəCj] as the driving force of schwa deletion. She did not identify the distributional pattern that decides whether schwa may be left out or not, though: schwa may be dropped when followed by any consonant plus yod, except if this consonant is the lateral. Has this generalisation gone unnoticed in the literature on [Cj] in French, then? 6 No: Dell (1973:262) incidentally notes that he will not address the question why schwa cannot be deleted when followed by a liquid and

V C V C V -C V C V C V C V -C V | | | | | | | | | | | | g i ʃ ə t j e o t ə l j e guichet-ier [ʃtj] hôtel-ier *[tlj]
By contrast under (4b), government cannot reach schwa because [lj] is not a good branching onset: the nucleus enclosed needs to be governed and therefore cannot be a governor itself. Thus schwa is pronounced.

Note that the issue of resyllabification does not arise in Strict CV since there are no branching constituents and there is no coda constituent (or post-rhymal complement): under (4a), the suffix-initial yod makes a branching onset with the root-final [t] by establishing an infrasegmental relation ( ).

Conclusion

The (Dell 1972:73;Tranel 1987:121;Klein 1992:23ff)? I pursue these issues elsewhere.

And there is also something to be understood about the whole [Cj] pattern in Quebec

French. Recall the parallel between [CəTR] and [CəCj] (Section 19.3.2): it holds indeed for the variety described in the article (unmarked Northern mainland or Parisian speakers), but Quebec

  Section 19.3.2: a number of speakers of Northern varieties cannot leave out schwa before TR). Selected speakers were asked to say for each word whether a pronunciation without schwa is possible, bizarre but possible, or impossible.

	) CəC-ier/-ière Results are unambiguous: items under (3b) where yod is preceded by the lateral were
	a. C1əC2j where C2 ≠ l
	b s gibecière overwhelmingly judged impossible by all speakers, while items under (3a) where a consonant 'gamebag' n t lunetier 'glasses maker'
	f t buffetier	'shopkeeper'	t panetier 'bread keeper'
	t cafetier different from the lateral precedes yod were overwhelmingly said to be either possible. Details 'café tenant' t panetière 'bread cupboard'
	t cafetière	'coffee maker'	t robinetier 'tap maker'
	l t giletier are as follows. 32 [CC+j] words (3a) and 15 [Cl+j] words (3b) were tested ([CCC+j] and [Cs+j] 'waistcoat maker' v chènevière 'hemp field'
	t giletière	'watch chain'	p t papetier 'paper maker'
	t molletière	'puttee'	r t charretier 'carter'
	t muletier	'donkey keeper'	t ferretier 'farrier hammer'
	t pelletier	'furrier'	t jarretière 'garter'
	t pelletiérine 'type of dewormer'	ʃ t guichetier 'counter clerk'
	t toletière	'oar reenforcer'	t n centenier 'centurion'
	m t cimetière	'cemetery'	v t buvetier 'barkeeper'
	n b canebière	'hemp field'	t louvetier 'wolf hunter'
	t bonnetier	'bonnet maker (masc.)'	t savetier	'cobbler'
	t bonnetière	'bonnet maker (fem.)'	z n dizenier	'commander of 10 men'
	t canetière	'silk worker'	t gazetier 'gazette owner'
	t chaînetier	'chain maker'	t noisetier 'hazel tree'
	t grainetier	'seed merchant'	s v sansevière 'sansevieria (bot.)'
	b. C1əC2j where C2 = l
	d l chandelier 'candle stick'	t l coutelier 'cutler'
	m l chamelier 'camel caravan driver'	l hôtelier	'hotel keeper'
	l sommelier 'wine waiter'	l râtelier	'rack'
	n l cannelier 'cinnamon tree'	l dentellière 'lacemaker'
	l tonnelier 'cooper'	s l chancelier 'chancellor'
	l tunnelier 'tunnel borer'	l pincelier 'brush box'
	p l chapelier 'hat-maker'	l vaisselier 'dresser'
	r l bourrelier 'harness-maker'	l ficelier	'trickery-using person'
	ʃ l bachelier 'bachelor'	l boisselier 'wooden objects merchant'
	l échelier	'ladder'	l muselière 'muzzle'
	l richelieu 'type of shoe'	l oiselier	'bird-seller'
	t l atelier	'workshop'	l roselière 'reed bed'
	l batelier	'ferryman'
	Grammaticality judgements by seven native speakers from mainland France that were
	born and brought up North of the Loire were collected. Only speakers were selected who can

leave out schwa in CəTR la secrétaire 'the secretary': this made sure that negative judgements of a particular item were not due to the impossibility for schwa to be left out before branching onsets (see

  yod ("La seconde question que nous laisserons de côté est le maintien de schwa lorsqu'il est suivi d'une liquide et d'un yod"). As far as I can see, he is the only author who has understood that schwa deletion in [C1əC2j] is blocked when C2 is the lateral. But since he says that he will not pursue this question, the pattern remains largely unillustrated in his book. It may thus be the case that the data in Section 19.4 are the first empirical demonstration of the generalisation at hand. Finally, note that Dell's statement also includes[rj] clusters: he says that schwa cannot drop in [C1əC2j] if C2 is a liquid, i.e. [r] or the lateral (note that [C1əC2]-ier/-ière under (3a) does not appear to produce any cases with C2 = [r]). It takes some argument to see that Dell is right (a piece of evidence are conditionals like donn-er-i-ez 'you (pl) would give' or huil-er-i-ez 'you (pl) would oil' where schwa cannot be dropped: *donn'rjez, *huil'rjez), but [rj] clusters represent an intricate pattern in French that cannot be examined here for lack of space.It was mentioned in Section 19.3.2 that Charette's 1991 analysis where[Cj] in[CəCj] is a branching onset in Quebec French is impossible in Standard GP since this would require resyllabification, which is prohibited in this theory: in cafet-ière for example, the lexically simplex onset containing the root-final consonant would have to become a complex onset in order to accommodate the yod coming from the suffix.

	19.6 Analysis based on branching onsets
	In Strict CV (Lowenstamm 1996; Scheer 2004), a branching onset identifies as a TR cluster
	belonging to two onsets whose solidarity is due to an infrasegmental relation (	in 4a) (Scheer
	2004: §64; Brun-Trigaud & Scheer 2010). The empty nucleus enclosed does not require
	government because it is silenced by this infrasegmental relation. Being therefore ungoverned,
	it is itself a good governor and governs the preceding schwa. This is shown under (4a) for [tj]
	in guichet-ier where schwa may be left unpronounced.
	(4) CəCj in Strict CV
	a. [Cj] branching onset: schwa governed b. [lj] no branching onset: schwa ungoverned
	Ggvt	Ggvt
	C

As we have seen, though, the pronunciation of schwa in [CəCj] depends on whether the following [Cj] cluster does (all [Cj] except [lj]) or does not

([lj]

) qualify as a branching onset.

  [CəCj] pattern is but a piece of the[Cj] puzzle in French, and the preceding pages have studied only a facet of that piece. Some other aspects and pieces are the following. Why is it that diérèse, i.e. expanding [j] into[ij], allows schwa to remain unpronounced in (pl) surround / bring to an end'), but impossible when originating in other suffixes (boucl-ier *[klj] 'shield')(Dell 1972:73)? How come that -i-ons/-i-ez can occur without [i] after [Tl] (vous boucl-i-ez [klj] 'you (pl) surround / bring to an end'), but not after [Tr] (vous plâtr-i-ez [trij] *[trj] 'you (pl) plaster')

conditionals (vous donn-er-i-ez [dɔn'rije] 'you (pl) would give'), but not in [Cəl]-ier (hôtel-ier *[ot'lije])? How come [Rlj] from /Rlj/ is possible (perl-ier [rlj] 'pearl, adj.', vous parl-iez 'you (pl) talked'), but not when [Rl] encloses a schwa /Rəlj/ (*bourr'lier 'saddler')? Why can schwa be absent in vous vous atteliez (att'liez [tlj]) 'you took care of', but not in atelier (*at'lier *[tlj])

'workshop'

(Klein 1992:39

)? Why is [Tlj] well-formed when created by -i-ons/-i-ez

(vous boucl-i-ez [klj]

'you