# Predicting Flowering Behavior and Exploring Its Genetic Determinism in an Apple Multi-family Population Based on Statistical Indices and Simplified Phenotyping

^{1}Laboratoire Jean Kuntzmann, Inria Mistis, Université Grenoble Alpes, Grenoble, France^{2}Virtual Plants Team, Inria and CIRAD, UMR AGAP, Montpellier, France^{3}Equipe Architecture et Fonctionnement des Espèces Fruitières, UMR AGAP, Institut National de la Recherche Agronomique, Montpellier, France^{4}Wageningen UR Plant Breeding, Wageningen University and Research, Wageningen, Netherlands^{5}Biometris, Wageningen University and Research, Wageningen, Netherlands^{6}Research & Technology Centre, Hendrix Genetics, Boxmeer, Netherlands

Irregular flowering over years is commonly observed in fruit trees. The early prediction of tree behavior is highly desirable in breeding programmes. This study aims at performing such predictions, combining simplified phenotyping and statistics methods. Sequences of vegetative vs. floral annual shoots (AS) were observed along axes in trees belonging to five apple related full-sib families. Sequences were analyzed using Markovian and linear mixed models including year and site effects. Indices of flowering irregularity, periodicity and synchronicity were estimated, at tree and axis scales. They were used to predict tree behavior and detect QTL with a Bayesian pedigree-based analysis, using an integrated genetic map containing 6,849 SNPs. The combination of a Biennial Bearing Index (BBI) with an autoregressive coefficient (γ_{g}) efficiently predicted and classified the genotype behaviors, despite few misclassifications. Four QTLs common to BBIs and γ_{g} and one for synchronicity were highlighted and revealed the complex genetic architecture of the traits. Irregularity resulted from high AS synchronism, whereas regularity resulted from either asynchronous locally alternating or continual regular AS flowering. A relevant and time-saving method, based on *a posteriori* sampling of axes and statistical indices is proposed, which is efficient to evaluate the tree breeding values for flowering regularity and could be transferred to other species.

## Introduction

Biennial bearing, defined as the irregular fruit or seed production over consecutive years, is a trait commonly observed in perennial crops (Monselise and Goldschmidt, 1982; Samach and Smith, 2013). In fruit trees, yield and fruit quality depend on bearing behavior, which is in turn strongly dependent on flowering intensity. However, floral induction may be inhibited by concurrent fruiting, leading to biennial bearing. Economically and environmentally sustainable techniques are therefore required for the management of biennial bearing in fruit production. An alternative strategy would be to select cultivars combining high fruit quality, long-term resistance to pests and diseases, tree architecture adapted to modern training systems and regular production. Breeding processes in trees are usually slower than in crops, because of the juvenile phase length and the long time required to assess the agronomic performances of pre-selected trees, especially bearing behavior (Laurens et al., 2012). Predicting bearing habit as soon as possible from the beginning of the genotype's production is thus of high interest. This strategy is reinforced by the existence of large differences among cultivars (Lauri et al., 1997, 2014) and the demonstration of genetic control of biennial bearing in an apple family derived from a cross between biennial and regular bearing parents (“Starkrimson^{®} Red Delicious” × “Granny Smith”; SG) (Guitton et al., 2012).

To characterize successive yields and bearing behavior, different approaches have been proposed (see Durand et al., 2013 for a review). The Biennial Bearing Index (BBI; see the list of abbreviations and notations in Table 1), which estimates the intensity of deviation in yields during successive years (Wilcox, 1944), has become the accepted standard to describe biennial bearing. It has been applied to yield (mass of fruit) at different scales: whole areas, individual trees or branches—on apple and other fruit tree species (Pearce and Doberšek-Urbanc, 1967; Reddy et al., 2003; Smith et al., 2004; Rosenstock et al., 2010; Guitton et al., 2012). However, the measure of the magnitude of irregular bearing by BBI is questionable, especially for trended series (Huff, 2001). A new methodology was thus introduced to characterize the bearing habit of trees as early as the first years of production, when the production is increasing. This methodology was based on a trend model on the yearly number of flowers, combined with a BBI-derived index and an index on correlations between residuals, denoted by γ (see Supplementary Material M1). An approximation of these indices based on within-tree sampling of successions of annual shoots (AS) along axes was considered. An entropy criterion was proposed to assess synchronicity of flowering in a given year, allowing a connection between axis- and tree-scale behaviors. However, the ability of axis-scale indices to predict genotype habits at tree scale was investigated on a single family (SG), and the annual shoot sequences were merely used to approximate the total number of flowering AS.

In the present study, we propose to extend the previous investigations by exploring new methods and indices based on the analysis of sequences of flowering shoots, and by performing a multi-family QTL detection to enlarge the genetic basis of biennial bearing variation in apple trees. We assumed that the analysis of entire sequences of successive AS, combined to flowering synchronicity in each year, would provide new insights on the genotype's behaviors. We thus proposed (i) to use not only the total number of flowering AS but also the vegetative ones and their succession; (ii) to derive new indices from this analysis at no additional measurement cost; (iii) to re-examine previous assumptions on the relation between alternation and regularity at tree and axis scales. Regarding genetics, we considered a larger germplasm to allow the comparison of alleles' performance in different genetic backgrounds (Pauly et al., 2012). We assumed that this will increase the number of segregating QTLs, detection power, accuracy of positions, and give more robust estimation of QTL effects (Bink et al., 2002; Liu et al., 2012). A Pedigree Based Analysis (PBA) was performed, using the concept of Identity By Descent (IBD) based on both pedigree and marker information (van de Weg et al., 2005; Luan et al., 2012). Our aim was to confirm previously found QTLs on the SG family and find new ones, tracing the original source of the favorable alleles and deepening our understanding of the genetic determinisms of biennial bearing in apple tree.

## Materials and Methods

### Plant Material

Five segregating families with known and related pedigree were used (Figure 1). The first family (Segura et al., 2006) is derived from a cross between a female parent (“Starkrimson^{®} Red Delicious”) having strong tendency to biennial bearing and a male parent (“Granny Smith”) prone to regular bearing. This family, hereafter referred as SG, is composed of 123 genotypes, each replicated twice in the same site, and among which 115 individuals were genotyped and phenotyped (Table 2). The second family, referred as XB (Celton et al., 2011), is from a cross between the hybrid X3263 (regular bearer) and “Belrène” (biennial bearer). It comprises 324 genotypes among which 50 were randomly selected for replication, resulting in a total of 374 trees that were all phenotyped, among which 58 were also genotyped. For both families, seedlings were grafted on a semi-dwarfing Pajam I rootstock and planted in 2004 and 2005, respectively, in a random experimental design at the Diascope INRA Montpellier experimental unit. All trees were grown under irrigated conditions with minimal training. In the SG family, branches along the trunk were removed below 50 cm in the first year and fruits were slightly thinned in the first 2 years of growth to avoid branch breaking. In XB family the trees were neither pruned nor the fruits thinned (Table 2).

**Figure 1. Genetic relationships between the five studied full-sib families (XB, HIVW, SG, N, and P; represented by black boxes) and their parents (represented by gray boxes) and founders or other members of the pedigree (represented by white boxes)**. Blue lines link the father to its progenies and red lines link the mother to its progenies. GoldenDel, “Golden Delicious,” ReiDuMans, “Reinette du Mans,” Wagenerap, “Wagenerapfel”; see text for family abbreviations and supporting information 2 for other abbreviations used in the pedigree. Reproduced from Allard et al. (2016) with the permission of Oxford University Press.

The three other families, called HIVW, N, and P, respectively, were chosen for their related pedigree. The HIVW family had a parent (X3263) common to the XB family and the other one (X3259) to the N family. Both N and P families had a common parent (X3305) and the three families derived from “Golden Delicious” with different parentage degrees (Figure 1). They are composed of 171, 42, and 45 individuals, respectively, each with a single replicate per genotype. They were planted at the INRA Angers experimental station, in 1992. The trees were trained in vertical axis with an annual manual thinning that was performed at the end of June and left one fruit per inflorescence. At both sites, pest and disease management was performed consistently with professional practices.

### Phenotyping

On the SG family, as described in Durand et al. (2013), successions of vegetative vs. floral AS were observed over consecutive years (based on the presence/absence of an inflorescence) along different types of axes: trunk, long and short axillary shoots (Table 2). Shoots were classified depending on their length. A distinction was also made between the long proleptic and sylleptic axillary shoots (see Segura et al., 2006). Flowering occurrence was observed along the trunk of each tree, as well as one long sylleptic and one long proleptic axillary shoot, both sampled on the first AS of the trunk (2004 AS). On each long axillary shoot, two short axillary shoots per AS were phenotyped the same way. Thus, 10 short axillary shoots of 5 to 1 years were recorded on long sylleptic, and 8 short axillary shoots of 4 to 1 years were recorded on long proleptic ones. The flowering pattern was described by recording the presence/absence of flowering event on AS (6 possible flowering occurrences on the trunk and long sylleptic shoots, 5 on the long proleptic axillary shoots). The data thus consisted in vegetative vs. floral AS in 6 to 1 year sequences, with 2,716 sequences in total, with a mean length of 3.0 (corresponding to 3 consecutive AS, or years, in average).

For the XB, HIVW, N, and P families, similar observations were performed on three long proleptic axillary shoots along which, four and three short AS were phenotyped along 2006 and 2007 AS, respectively. The numbers of sequences were 7,757 in XB (with mean length of 5.3), 1,511 in HIVW (mean length of 6.0), 442 in P (mean length of 6.2), and 905 in N (mean length of 6.4; see distributions in Supplementary Figure S1).

### Statistical Modeling of AS Fate Sequences

Our approach is based on the classical BBI and on indices defined in Durand et al. (2013): BBI-derived indices (denoted BBI_norm and BBI_res_norm), autoregressive coefficient γ_{g} and entropy. They are based on counts of flowering AS at axis and whole-tree scales, whenever possible. The description of the trend and autoregressive models, BBI-derived indices and the statistical methodology for classification of genotype habit can be found in Supplementary Material M1. Compared to the original indices, we added a fixed “site” effect, Montpellier (M) or Angers (A), in the trend and autoregressive models, whenever possible.

At the axis scale, sequences of AS fates are denoted (*F*_{g, r, π, t, ℓ})_{ℓ≥0} with (*F*_{g, r, π, t, ℓ} = 0) denoting the absence and (*F*_{g, r, π, t, ℓ} = 1) the presence of flower for replication *r* of genotype *g* at site π, year *t*, and location (or AS) ℓ in the axis. The indices at axis scale, denoted by BBI_ax, BBI_norm_ax, BBI_res_norm_ax, and γ^{ax}, were computed as those defined at tree scale but using the total yearly counts of flowering AS in axes sampled within each tree replicate.

These counts were also used to compute two entropies, denoted respectively ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$(entropy based on frequencies) and ${\overline{\begin{array}{c}Ent\end{array}}}_{glmm,g}$ (entropy based on a generalized linear mixed model-GLMM). These two indices are based on the assumption of independent Bernoulli distributions for the successive AS fates *F*_{g, r, π, t, ℓ}. This assumption led us to ignore dependencies and patterns of alternation that could be inferred from sequences of AS fates along axes.

It can be assumed that models and indices taking explicitly into account the succession of AS fates would yield better predictions of genotype habit and provide new insights on the relationship between alternation at axis and whole tree scales. This is why we modeled the sequence of AS fates (*F*_{g, r, π, t, ℓ})_{ℓ≥0}, in this new study, using high-order Markov chains (see Costes and Guédon, 2012, for discussion of the shortcoming of using first-order rather than high-order Markov chains). In such models, the variable at time *t* depends on the *M* past variables where *M* is the order or memory length. The values of these *M* past variables are referred to as the memory *m* of the model at time *t*. For example, memory “10” means that flowering occurred at time *t*–2 but not at time *t*–1. Markov chains with different orders were estimated on the basis of every non-overlapping sequence extracted from the trees. A second-order Markov chain was chosen by a model selection procedure based on the Bayesian Information Criterion (BIC, see Kass and Raftery, 1995). This choice of order means that knowledge of the presence/absence of flowers at years *t*–1 and *t*–2 is necessary for prediction of flowering at year *t*. Thus, the set of memories was {00, 10, 01, 11}.

Assuming that alternation is partly genetic, some interactions between year *t*, memory *m* and genotype *g* should have an effect on flowering. To model these interactions with binary observations *F*_{g, r, π, t, m, ℓ}, approaches based on GLMMs are relevant (Molenberghs and Verbeke, 2005). The following Markovian GLMM was considered:

where λ is a fixed intercept, ρ_{π} is the fixed effect of site π (with Angers site as reference ρ_{A} = 0), μ_{m} is the fixed effect of memory *m* (with reference μ_{00} = 0), φ_{t} is the fixed effect of year *t* (with reference φ_{2006} = 0) treated as a qualitative variable, variables θ_{g, m} are independent random interactions between genotype *g* and memory *m* with common variance ${\tau}_{\theta}^{2}$, variables η_{g, t} are independent random interactions between genotype *g* and year *t* with common variance ${\tau}_{\eta}^{2}$, and variables ζ_{g, r} are independent replication-specific random effects with common variance ${\tau}_{\zeta}^{2}$. All random effects were assumed to be mutually independent and Gaussian. This model consists of a high-order Markov chain for process (*F*_{g, r, π, t, m, ℓ})_{ℓ≥0} where the transitions are treated as GLMMs, so as to introduce fixed and random effects in modeling binary outcomes. Parameter estimation was by restricted maximum likelihood. For a better interpretation of the model and obtain new indices, it is useful to estimate the value of the random effects. This is achieved by their Best Linear Unbiased Predictors (BLUPs), which are the conditional means and also the most probable values of the random effects. These were computed by *glmer* and *ranef* functions of package *lme4* (Bates et al., 2011). The BLUPs of θ_{g, m} were used to discriminate genotypes on their low vs. high probability of AS bearing flowers at year *t* given they had memory *m*. Similarly the BLUPs of η_{g, t} were used to discriminate genotypes on their low vs. high probabilities of bearing flowers at year *t*.

### Predicting Tree Fruiting Behavior from Axis-Scale Indices

One goal was the prediction of tree flowering behavior with respect to three classes: regular, irregular or biennial. We assumed that the true class of each genotype could be deduced from the tree-scale indices, for the SG family, using a clustering method developed in Durand et al. (2013). The dependencies between class and tree-scale indices are represented in Figure 2. The classification of the genotypes of all families from axis-scale indices was achieved using as predictors: the BLUPs of the θ_{g, m} random effects, BBI_res_norm_ax, γ^{ax}, ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$and ${\overline{\begin{array}{c}Ent\end{array}}}_{glmm,g}$ (which were available for every genotype, as opposed to whole tree-scale predictors). Note however that BLUPs may not be defined for every parameter and genotype. For example, if memory 11 did not occur in the axes of some genotype *g*, θ_{g, 11} cannot be defined, resulting into missing predictors. Moreover, some predictors could be highly correlated and redundant. To handle both issues of predictor absence and redundancy, a principal component analysis (PCA) for partially missing data was used, with the R package missMDA (Josse and Husson, 2012). Classification was performed using neural networks (NNs, see Supplementary Material M2). The number of principal components (PCs) and the NN regularization parameter were determined by out-of-sample validation.

**Figure 2. Representation of the clusters for the genotypes in the SG family, as a function of the tree-scale indices BBI_res_norm and γ _{g} (auto.cov)**. Cluster 1 can be interpreted as regular bearing genotypes, cluster 2 as biennial bearing genotypes, and cluster 3 as irregular bearing genotypes. Reproduced from Durand et al. (2013) with permission from Oxford University Press, Copyright 2013.

Although classification can be a relevant method to assess regularity at tree scale, its prediction is a rough summary (through three classes only) of the flowering behavior. In other contexts, such as genotype ranking, quantitative assessment of the bearing behavior may be required. This could be achieved through the tree-scale indices BBI_res_norm and γ_{g} (when measured), since they provide a more accurate description of this behavior. Since tree-scale indices were known for SG only, approximation of both indices was performed by regression, from axis-scale indices θ_{g, m}, BBI_res_norm_ax, γ^{ax}, ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$and ${\overline{\begin{array}{c}Ent\end{array}}}_{glmm,g}$. NNs were used as nonlinear regression functions (instead of nonlinear classifiers in the case of classification), also using missing data PCA. The NN parameters were estimated by least squares minimization. Since the optimal numbers of PCs to be used in classification and regression NNs may be different, they were both chosen independently by out-of-sample validation. The approximated values of BBI_res_norm and γ_{g} are referred to as BBI_res_norm_pred and γ^{pred}, respectively. These two values are necessarily linearly dependent due to the nature of NNs.

Classification and regression were also performed in Durand et al. (2013) but the models did not consider ${\overline{\begin{array}{c}Ent\end{array}}}_{glmm,g}$ nor θ_{g, m}. Therefore, we assessed the gain of using these new indices as predictors. Since the tree-scale indices were known for SG only, the classification and regression errors were assessed on this family.

### Genetic Map and QTL Mapping

The five full-sib families and their progenitors were genotyped with the Infinium^{®} 20K SNP array (Bianco et al., 2014), according to Chagné et al. (2012) and Antanaviciute et al. (2012). A genetic map composed of 7,100 SNPs has been integrated over 27 full-sib families using the same approach as in Di Pierro et al. (2016), and was used for QTL mapping (van de Weg et al., unpublished). 6,849 SNPs were used after careful checking of their robustness (van de Weg et al., 2013; Di Guardo et al., 2015), consistency and recombination pattern on the 5 families and the pedigree members (Allard et al., 2016). The quality of the map was achieved by intense data curation and by using graphical genotyping to avoid double recombinations along with the use of multiple families to create an integrated genetic map that reduced cases of false marker order. The high quality of our current map is underlined by the low number of SNPs that are in discordant order (71 SNPs, 3.2%), the small size of the genetic segments in which these discordant orders occurred (usually 0.5 cM, data not shown), and the similar small size of both genetic maps. Then sets of single SNPs were integrated into haploblocks, corresponding to successive 1 cM segments. Haplotypes were composed using the software FlexQTL™ (www.flexqtl.nl) and PediHaplotyper (Voorrips et al., 2016).

Indices used in QTL analysis were BBI_res_norm_ax, γ^{ax}, BBI_res_norm_pred, γ^{pred}, ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$, ${\overline{\begin{array}{c}Ent\end{array}}}_{glmm,g}$, the BLUPs for genotypes × memory interactions (θ_{g, 00}, θ_{g, 01}, θ_{g, 10}, θ_{g, 11}) and genotype × year interactions (η_{g, 2006}, η_{g, 2007}, η_{g, 2008}, η_{g, 2009}, η_{g, 2010}, η_{g, 2011}, η_{g, 2012}). Among all variables, QTLs were detected using a linear model that comprised an intercept μ, the regression coefficients *a* on the QTL covariates, and a residual *e*, as:

where *W* is the design matrix for the QTL effects. A bi-allelic model is assigned to a QTL with alleles denoted by Q and q, with only additive effects and values of [QQ, Qq, qq] equal to [1, 0, −1]. As QTL genotypes of individuals are a priori unknown, modeling is based on independent assignment of alleles to founders and segregation indicators to trace transmission from parents to offspring (Bink, 2002). Uniform priors were assigned to μ and λ (vector giving QTL position), and normal priors to the vectors *a* and *e* in (2), i.e., *a*~$N(0,{\sigma}_{a}^{2}I)$ and *e*~$N(0,{\sigma}_{e}^{2}I)$. ${\sigma}_{a}^{2}$ and ${\sigma}_{e}^{2}$ are the per-QTL explained variance and the residual variance, with priors being inverse Gamma distributions (Bink et al., 2008). The number of QTLs was assigned a Poisson prior. Results for a prior mean of 5 are reported only. Other values yielded similar results and inferences (data not shown). Samples from the joint posterior distribution $f(\mu ,a,\lambda ,{\sigma}_{a}^{2},{\sigma}_{e}^{2}|y)$ of the model parameters were obtained by Markov chain Monte Carlo (MCMC) simulation in FlexQTL™ (Bink et al., 2008, 2014).

The MCMC algorithms and details on the monitoring of Monte Carlo accuracy and length of the simulation chains can be found in Bink et al. (2008) and Allard et al. (2016), respectively. The number of QTLs was inferred from a pairwise comparison of models differing by one QTL, and considering twice the natural logarithm of the Bayes Factors (Kass and Raftery, 1995), denoted 2*lnBF. Values >2, 5, and 10 indicate positive, strong, and decisive evidence for the presence of a QTL, respectively. QTL positions were based on posterior QTL intensities, and QTL contributions on the posterior mean estimates of the QTL effects. Posterior probabilities of QTL genotypes were also estimated (Bink et al., 2014).

When several QTLs were detected for a variable, the interactions between QTLs were tested by linear models with haplotypes located at the peak of the QTLs. Model selection was achieved with a backward method based on AIC (Burnham and Anderson, 2002).

## Results

### Modeling as Sequence

A Markovian GLMM was estimated merging the five families (1) in Section Statistical Modeling of AS Fate Sequences). It had a BIC value of 44,300. It was compared with the models (i) without any random effect (i.e., with ${\tau}_{\theta}^{2}={\tau}_{\eta}^{2}={\tau}_{\zeta}^{2}=0$, BIC = 52,944), (ii) containing “genotype” and “replication” random effects only (no interactions with year or memory, BIC = 50,723), (iii) without “replication” random effects (BIC = 44,394) and (iv) without “genotype” random effects and their interactions (BIC = 50,780). The parameter estimates are in Table 3.

**Table 3. Estimates of fixed effects and variances (with the p-values of the tests of the null hypothesis “n = 0” against the alternative “n ≠ 0,” for parameters n associated with fixed effects) of the mixed model estimating the probability of flowering at axis scale in five apple tree families**.

The BIC values of these models showed that all fixed (site, memory, year) and the random replication effects were significant, consistently with the associated *p*-values (Table 3). The replication effect, although included in the best model, induced less variability in flowering than memory and year effects (the latter having the largest variance). The trees in Montpellier had lower flowering probability than those in Angers. Flowering AS were more frequent after a vegetative AS preceded by a flowering AS, than directly after a flowering AS (since μ_{10} was higher than μ_{01} and μ_{11}), showing frequent biennial alternation in flowering at axis scale. The probability of flowering was the highest in 2008, whereas it was particularly low in 2010 and 2012, whatever the site.

Empirical standard deviations were computed for the BLUPs of random effects θ_{g, m} and η_{g, t} to estimate their specific variability for each memory *m* or year *t* (Table 4). The random interactions θ_{g, 10} (between genotypes and memory “flowering AS followed by a vegetative AS”) had the lowest genetic variability. In contrast, the random interactions θ_{g, 01} between genotypes and memory “vegetative AS followed by a flowering AS” had noticeably higher variability. The genetic variability of the random interactions η_{g, t} increased with years, showing that the genetic differences became larger with tree age.

**Table 4. Empirical standard deviations of random effects for each kind of interaction of the mixed model estimating the probability of flowering at axis scale in five apple tree families**.

### Prediction of Tree-Scale Indices Using Linear Regression of Axis-Scale Indices

Trend models of the yearly numbers of flowers at axis scale were used to compute BBI_res_norm_ax and γ^{ax} [models (A) and (I) in Supplementary Material M1]. However, the site effect was included in (I) only, since its inclusion in (A) induced non-identifiability issues.

BBI_res_norm and γ_{g} were regressed with three or four principal components (PCs) using NNs, and the best cross-validated correlations were obtained with three PCs. The optimal correlation between BBI_res_norm and its prediction BBI_res_norm_pred was 0.71 when using BBI_res_norm_ax, γ^{ax} and ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$ as predictors, and 0.72 when using the PCs. The optimal correlation between γ_{g} and γ^{pred} was 0.60 using BBI_res_norm_ax, γ^{ax} and ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$, and 0.64 using PCs. Using ${\overline{\begin{array}{c}Ent\end{array}}}_{glmm,g}$ instead of (or in addition to) ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$ did not improve correlations. Thus, adding information from Markovian GLMMs did not significantly improve the prediction of the tree-scale indices.

### Classification of the Genotypes with Respect to Tree-Scale Bearing Habit

Classification of the genotypes of SG, which may be interpreted as the expected error rate on other families, yielded a 39% cross-validated error rate when using BBI_res_norm_ax, γ^{ax} and ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$, and 35% with PCs. Genotypes with unknown values of BBI_res_norm_ax, γ^{ax} or ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$were excluded, thus keeping 115 genotypes in SG. As previously, the best prediction was achieved with three PCs. Both ${\overline{\begin{array}{c}Ent\end{array}}}_{glmm,g}$ and ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$ were included in the PCA, in addition to the other axis-scale indices. Although limited, the improvement of the error rate, was significant at level 0.7% on 50 random test samples (Student's *t*-test).

To predict the unknown bearing habits at tree-scale of the genotypes of XB, HIVW, N, and P, the optimal NN model on SG was re-estimated on the whole data set, using genotypes with known classes (i.e., 122 genotypes in SG) for learning the mapping between local indices and classes. Confusion between classes concerned irregular genotypes that could hardly be discriminated from regular and biennial bearing genotypes (Table 5). This comes from irregular genotypes having intermediate values of their indices, between those for the regular and biennial bearing genotypes. As a result, 15 regular and 9 biennial genotypes were classified as irregular, and 9 irregular genotypes were classified as regular. In contrast, one misclassification only occurred between regular and biennial bearing genotypes, highlighting that discrimination between both behaviors is easy.

**Table 5. Contingency table for the number of genotypes of each possible true class (corresponding to observations on SG family) assigned to each possible predicted class by NN on local indices**.

The three classes were discriminated by analyses of variance (ANOVAs) performed on the axis-scale indices (see Supplementary Table T1). The lowest *p*-values (<1e-8) were obtained with θ_{g, 01}, ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$, BBI_res_norm_ax and γ^{ax}, which pointed out especially contrasted values of these indices among the three classes. These four indices thus highlight a higher potential than the other ones to discriminate between classes of bearing behavior.

The model yielded the following predictions:

• SG: 29 (24%) regular, 29 (24%) biennial bearing, 64 (52%) irregular genotypes—with 33% error rate on the learning sample (37% if using BBI_res_norm_ax, γ^{ax} and ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$ only)

• XB: 48 (17%) regular, 69 (25%) biennial bearing, 158 (57%) irregular genotypes

• HIVW: 38 (22%) regular, 23 (14%) biennial bearing, 109 (64%) irregular genotypes

• N: 31 (31%) regular, 9 (9%) biennial bearing, 60 (60%) irregular genotypes

• P: 10 (19%) regular, 4 (8%) biennial bearing, 39 (74%) irregular genotypes.

### Consistency of the Indices between Tree and as Scales

Correlation matrices were computed between the tree and axis-scale indices for SG (Table 6) and between indices at axis-scale only for the five combined populations (Table 7). The strongest correlations between tree-scale and axis-scale indices were the correlations between BBI_res_norm_ax, γ^{ax} and θ_{g, 01} (which correlation is negative with BBI_res_norm_ax and positive with γ^{ax}). Moreover, θ_{g, 11} and ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$ were moderately correlated with BBI_res_norm (negatively) and with γ^{ax} (positively), while ${\overline{\begin{array}{c}Ent\end{array}}}_{glmm,g}$ was poorly correlated and θ_{g, 10} uncorrelated (at level 0.05) with the tree-scale indices.

**Table 6. Correlation coefficients between indices at whole tree scale and indices and BLUPs at axis scale, with 95% confidence intervals, in the SG family**.

**Table 7. Correlation coefficients between the genotype × memory interactions (θ _{g,m}) BLUP indices and the other indices and BLUPs at axis scale, with 95% confidence intervals, combining the five apple tree families**.

### QTL Mapping

QTLs were detected for all indices except for θ_{g, 10}, η_{g, 2007}, η_{g, 2009}, η_{g, 2010}, η_{g, 2011}, and η_{g, 2012}. Two major QTLs were detected with a strong evidence (2*lnBF ≥ 5) on LG4 and LG, for the two BBI-derived indices (BBI_res_norm_ax, BBI_res_norm_pred) (Table 8, Figures 3A,B, Figure S2 in Supplementary Material for the trace plot of QTL across iterations). The QTL detected on LG4 explained 11.5, and 13.3% of variance, respectively. The QTL on LG5 explained 6.9, and 8.3% of the variance of each index, respectively (Table 8). Two other QTLs were detected on LG8 and LG10 but had a strong evidence for BBI_res_norm_pred only. They explained 11.7 and 10% of variance of this index, respectively.

**Table 8. Parameters associated with the QTL detected for BBI derived indices, autoregressive coefficients and entropy**.

**Figure 3. Posterior probability of QTL position along genome, the beginning and the end of the chromosomes are represented by vertical dashed lines**. The variables displayed are **(A)** BBI_res_norm_ax, **(B)** BBI_res_norm_pred, **(C)** γ^{ax}, **(D)** γ^{pred}, **(E)** ${\overline{Ent}}_{g}$, **(F)** ${\overline{Ent}}_{glmm,g}$. See text for abbreviation meaning.

Two QTLs were detected on LG8 and LG10 but had a strong evidence for BBI_res_norm_pred only. They explained 11.7 and 10% of variance of this index, respectively.

For the autoregressive coefficients, the same regions along the genome were detected. For γ^{ax}, three QTLs were detected on LG4, LG 5 and LG 10 (Figure 3C) but none with a strong evidence (Table 8). The QTLs on LG5 and LG10 colocalized with that of the BBI indexes (Figure 3). Four QTLs were mapped for γ^{pred} on LG4, LG5, LG8 and LG10 (Figure 3D) and colocalized with BBI_res_norm_pred. QTLs on LG4 and LG8 explained 12%, and those on LG5 and LG10 explained 10% of the variance of γ^{pred} (Table 8). Five QTLs were detected for ${\overline{\begin{array}{c}Ent\end{array}}}_{glmm,g}$ (Figure 3F), among which only that on LG9 had a strong evidence and explained 20% of the variance.

The QTL detected for θ_{g, 00}, θ_{g, 01} and_{θg, 11} (Table T2; Figures S3, S4 in Supplementary Material) on LG10 colocalized with that detected for BBIs. It explained 10, 11.8, and 10.9% of the variance, respectively, and had a strong evidence for θ_{g, 00} and θ_{g, 01}. Two QTLs with positive evidence were mapped on LG12 and LG6, for η_{g, 2008} and η_{g, 2006}, respectively (Table T2; Figures S3, S4 in Supplementary Material).

No interaction between QTLs could be identified for BBI_res_norm_ax and γ^{ax}, whereas interactions were identified for BBI_res_norm_pred between QTLs on LG4 and LG8, and for ${\overline{\begin{array}{c}Ent\end{array}}}_{glmm,g}$ between QTLs on LG7, LG9, and LG17, with all 2-way and 3-way-interactions being significant (results not shown).

### Genotype Estimation at Main QTLs

Genotype estimation at QTLs brought two types of information. Firstly, the allelic classes of parents allowed identifying in which family QTLs segregated: if one parent of a family was estimated to be heterozygous at a QTL, this family segregated for this QTL. Secondly, genotype estimation allowed identifying parents and founders bearing favorable alleles, specific to the considered variable. Hereafter, only QTLs with strong evidence were commented (Table 9).

**Table 9. Parent genotype estimation for BBI_res_norm_ax at each QTL: qq for homozygous with low value favorable allele, QQ for homozygous with high values unfavorable allele, and Qq for heterozygous**.

The estimated genotypes were identical for all QTLs that colocalized for the normalized BBI indexes, except for that on LG5 for BB_res_norm_pred (Figure 3). The LG4 QTL segregated in SG only and was heterozygous for “Red Delicious.” The other parents were estimated to be homozygous for the low value allele. As low BBI values indicate a regular bearing habit, “Red Delicious” is likely to transmit the favorable allele to half of its progenies, whereas the other parents are supposed to transmit favorable alleles only.

The QTL on LG5 segregated in all families except XB with parents X-3305, X-3259, and “Red Delicious,” which is estimated to be heterozygous. However, for BBI_res_norm_pred, X-3259 was estimated to be homozygous (Table 9), this questioning the presence of favorable allele at that position. “Rubinette” was estimated to be homozygous for the low value allele therefore supposed to transmit only favorable allele. “Granny Smith” was estimated to be homozygous for the high value allele therefore supposed to transmit unfavorable alleles only.

Three parents, “Rubinette,” “Granny Smith,” and “Belrene,” were estimated to be heterozygous at the LG10 QTL, consistently with the contribution of several families to this QTL (Figure 3). Unstable results were found for the genotype estimation of X-3263, depending on the index used.

“Granny Smith” was estimated to be heterozygous on the LG9 QTL for ${\overline{\begin{array}{c}Ent\end{array}}}_{glmm,g}$, consistently with the segregation of this QTL in SG only. The other parents were estimated to be homozygous for the high value allele (Figure 3). Since individuals with the most regular fruiting behavior had a high value of entropy, these homozygous parents are supposed to transmit favorable alleles.

The genotype estimates over all QTLs for a given trait provide information on the value of the parents and founders as genitors. Taking BBI_res_norm_ax as an example, antagonisms between QTLs were revealed: the parents estimated to be homozygous for the favorable low value allele on LG4 were estimated to be almost systematically homozygous for the unfavorable allele on LG7 (Tables 8, 9). The count of favorable alleles over the six QTLs revealed that the founder “Coop17” and the parent X-3259 had the highest values (12 and 9 favorable alleles, respectively). In spite of QTLs for which no genotype could be estimated, “Rubinette” and “Golden Delicious” appeared interesting parents with homozygous estimated genotypes with the favorable allele at four QTLs over six (Table 9).

## Discussion

### Efficiency of Indices Derived from Whole as Sequence Analyses at Tree Scale

In this study, we estimated several indices to capture the bearing behavior of a large set of genotypes. As previously underlined (Durand et al., 2013), BBI_res_norm and γ_{g} are negatively correlated and complement one another: BBI_res_norm distinguishes between regular individuals, with low values, and irregular and biennial bearing individuals, with high values. As a complement, γ_{g} distinguishes biennial bearing individuals, with negatives values from regular and irregular individuals. The new indices θ_{g, m} and their correlations with γ_{g} and BBI_res_norm provide complementary information: the regular genotypes (lowest values of BBI_res_norm) exhibit AS with a flowering probability above average at year *t* after flowering at year *t*–1 (highest θ_{g, 01} and θ_{g, 11}, associated with memories 01 and 11). In contrast, θ_{g, 00} and θ_{g, 10} are less discriminant, because the probability to flower after vegetative events is always high. The positive correlation of entropy ${\overline{\begin{array}{c}Ent\end{array}}}_{g}$ with γ_{g} and its negative correlation with BBI_res_norm show that the genotypes with the highest synchronism (lowest value of entropy) are mostly biennial bearers (high values of BBI_res_norm, low values of γ_{g}).

The correlations between tree- and axis-scale descriptors suggest that biennial bearing at tree scale results from the conjunction of two phenomena: synchronism in flowering between AS in a given year and biennial alternation at AS scale between consecutive years. On the contrary, regularity at tree scale results from either asynchronous locally alternating flowering or regular flowering at AS scale. Irregular genotypes exhibit intermediate values for every descriptor, suggesting that these genotypes are characterized by partial biennial alternation at AS scale or strong biennial alternation with partial synchronism. However, more complex within-tree organization of synchronisms could exist (Couranjou, 1983), that have not been investigated herein. The regular genotypes exhibiting synchronized and regular flowering AS or desynchronized flowering AS will require further investigations regarding fruit set and quality. Indeed, high flowering rate is usually associated with poor fruit set due to environmental (Tustin et al., 2012) or genetic co-variation (Celton et al., 2014). Tree management includes a number of practices to reduce crop load, such as thinning (manual or chemical) or the manual removal of fruiting spurs (Lauri et al., 2007; Breen et al., 2016). We can thus suspect that selecting genotypes with regular desynchronized axes could be an appropriate strategy for avoiding poor fruit set while reducing thinning or manipulations costs. In conclusion, three indices can be considered as key and complementary descriptors of the bearing behavior of genotypes at either tree or axis scales: BBI_res_norm, γ_{g} and entropy. The first two are sufficient to classify the genotypes into regular, irregular and biennial classes. However, entropy allows this diagnostic to be refined by providing information on the within-tree strategy of regular genotypes, with potential consequences on tree management and breeding goals.

More insight on the bearing behavior is also gained by introducing site and year effects and analyzing the genotype × year interactions (η_{g, t} indices). The lower flowering probability at Montpellier than Angers may result from the absence of thinning practices on XB family, which may have hampered tree flowering capacity over years (Dennis and Neilsen, 1999). However, thinning was performed quite lately in Angers, due to a large dispersion of phenological stages among genotypes (Allard et al., 2016). This practice likely had a relatively low impact on floral induction which is assumed to occur mid-June in apical meristems of spurs (Hanke et al., 2007) and therefore on alternation. Even though the mean probability of flowering of all genotypes per site highlighted phase opposition in the last three years (2010–2012, Supplementary Figure S5, left), no clear characterization of years as being “ON” or “OFF” could be made on the mean values per family (Supplementary Figure S5, right). Thus, no climatic year could be considered as “ON” or “OFF.” Even though critical climatic conditions such as frost (Nagy et al., 2010) or crop load management (Girona et al., 2010) can synchronize trees in a given year and site, different genotypes can be in phase opposition for flowering in a given year, in apple (Durand et al., 2013) as in olive tree (Ben Sadok et al., 2013).

The new indices at axis scale were expected to improve prediction of bearing habits at tree scale. Even though the predicted BBI_res_norm and γ_{g} appeared robust based on QTL detection, only a 4% improvement in predictions on test samples was obtained compared with Durand et al. (2013) and the classification error was still of 35% (even though that on SG family was reduced from 37 to 33%). The misclassification mainly concerned the irregular genotypes, whereas the regular and biennial behaviors could be predicted with good accuracy. The misclassification of regular genotypes considered as irregular (15 over 36, see Table 5) could lead to discard them during the selection process. However, this type of error is less problematic than the reverse (selecting irregular genotype that would be misclassified as regular) especially if we consider the drastic reduction in the number of individual selected in the early stages of breeding process. Therefore, the simplified phenotyping strategy that consists in sampling axes within the tree structure with an *a posteriori* observation appears to be relevant. Indeed, this is a time-saving strategy for phenotyping that can be combined to computation of indices for rejecting biennial and irregular genotypes during the assessment of agronomic performance of pre-selected genotypes in breeding programmes. It could also enable further phenotyping of germplasm toward the implementation of DNA-informed breeding approaches by further enlarging the number of founders and breeding parents for which QTL-genotypes are known, or by implementing genomic selection. Both applications are likely to accelerate the breeding progress and overcome long generation intervals and extensive phenotyping in outbred fruit tree crops (Kumar et al., 2013; Muranty et al., 2015).

Moreover, this strategy could be directly used on any species for which retrospective phenotyping of flowering is possible at AS scale. This is the case for species with terminal flowering such as pear, walnut, avocado, mango, litchi etc. in which flowering events can be easily identified. For such species, the methodology proposed, including prediction of flowering behavior at tree scale from *a posteriori* observations and computation of indices would be transposable. For other species, BBI_res_norm and γ_{g} index could be computed based on counting the total number of inflorescences measured on several successive years. Even though more time-demanding than retrospective observations, such counts may be facilitated and automatized by new technologies based on imagery (Aggelopoulou et al., 2010; Gongal et al., 2016).

### Genetic Determinisms of Bearing Behavior in a Multi-family Population

Five major QTLs were yielded, four for BBIs and auto-correlation coefficient (γ_{g}), on LG4, 5, 8, and 10, and one for entropy on LG9, which were partially common with previous studies on the SG family, such as the QTL on LG4 previously detected for BBI (Guitton et al., 2012). Also, a QTL on LG8 was found for BBI at both tree and axis scales in Durand et al. (2013) and for BBI, yield and number of flowers per inflorescence in Guitton et al. (2012). This zone, located at 8–23 cM on LG8, partially overlapped with QTLs detected in SG family for a descriptor of tree vegetation density (Virlet et al., 2015), for traits linked to bud break (Celton et al., 2011; Allard et al., 2016), and traits involved in gas exchange and xylem conductance (Segura et al., 2008; Regnard et al., 2009; Lauri et al., 2011). The QTL on LG10 located between 55 and 78 cM for BBIs, γ^{ax}, θ_{g,00} and θ_{g,01} co-localized with those detected for BBI, precocity and number of seeds per inflorescence in SG (Guitton et al., 2012) and for the percentage of bourses with one fruit on short axes in XB family (Celton et al., 2014). The QTLs detected on LG1 and LG14 by Durand et al. (2013) were confirmed, but in year specific interaction only, whereas the QTL on LG11 could not be confirmed. Actually, re-analyzing the same dataset led us to found an inappropriate account for missing flowering AS that had led to a false QTL detection. This was corrected in the present study.

As previously suggested (Bink, 2002; Liu et al., 2012), a higher power of detection was obtained in a multi-family context, which brought a higher number of segregating regions, alleles and individuals. The QTL on LG9 for entropy appeared as a new zone of importance, as it co-localized with a major QTL detected for the timing of vegetative and flowering bud break (Dyk et al., 2010; Celton et al., 2011; Allard et al., 2016). New QTLs were also detected for BBIs on LG5 and LG7, consistently with “Starkrimson^{®} Red Delicious” and “Granny Smith” not being heterozygous for these QTLs. The QTL on LG5, located from 9 to 24 cM, co-localized with QTLs previously detected for variables linked to the tree fruiting capacity (number of fruits and fruit biomass) under soil water restriction (Virlet et al., 2015). Moreover, LG5 is homologous of LG10 (Velasco et al., 2010; Bushakra et al., 2012), also involved in the fruiting capacity of the trees. As LG5 and LG10 are full-length homologs which orientation is defined upside-down, these two QTLs may have a common underlying mechanism. Altogether these co-localizations suggest that both tree development (LG8) and fruiting capacity (LG5 and LG10) may contribute to the genetic variation of biennial bearing behavior in apple tree. They reinforce previous assumptions regarding the combined effects of both competition among organs for nutrients and hormonal signals on biennial bearing (e.g., Chan and Cain, 1967; Dennis and Neilsen, 1999). To further decipher the putative role of fruits and carbon economy on the inhibition of floral induction, the tools defined herein could be wisely used to classify genotypes before investigating their physiological behaviors.

The lack of QTL interaction detected for BBI_res_norm_ax and γ^{ax} suggests that QTL contribution to the genetic variance was properly estimated for these variables. However, the QTL interactions detected for BB_res_norm_pred and ${\overline{\begin{array}{c}Ent\end{array}}}_{glmm,g}$ suggest that the QTL contribution might be underestimated. Taking into account epistasis might provide a better understanding of the genetic architecture of these traits. The complex genetic architecture of the studied traits and the alleles present on the different QTLs may lead to different degrees of alternation. Further characterization of allelic variations will be necessary for analyzing their relative contribution to the tree phenotype.

Finally, no co-localization was found with QTLs detected for architectural traits measured in the first years of tree development (Segura et al., 2007, 2009). Even though qualitative notations of architectural traits collected on young trees and their linear combinations could lead to an early diagnostic on biennial bearing (Lauri et al., 2014), no correlations or co-localization could be found in SG family, which was the only one studied for early tree architectural development among the five families considered here. In future work, including the type of bourse shoot within successive floral AS could improve the characterization of genetic variations and their relationships with architectural factors.

### Potential Use in Breeding of Genitors or Founders

In a breeding perspective, three descriptors should be combined to ensure regular bearing behavior (i.e., BBI_res_norm_pred, γ^{pred} and entropy). Five major QTLs on LG 4, 5, 8, 9, and 10 should be checked and alleles adequately combined in new released materials. Considering the strong evidence of QTLs for these three descriptors, appropriate phenotypes could be targeted with low BBI_res_norm_pred and medium or high γ^{pred} at the tree scale with high entropy values avec the axis scale. By contrast, trees with low values of BBI_res_norm_pred, medium or high γ^{pred} values and high entropy values at the tree scale could not be observed in the studied populations. As underlined by Samach and Smith (2013), the evolutionary advantage of masting (i.e., synchronicity of flowering at tree and population levels) remains questionable. Our results suggest that flowering synchronicity at the tree level could not be associated with regularity probably because it would lead to over-cropping and major drawbacks in an agronomic context. Knowing if flowering desynchronization has been selected during the apple domestication remains an open question.

This study revealed a complex genetic architecture of flowering habit in apple. The overview of all loci involved in trait variation led us to assess promising individuals and progenitors. X-3259 appeared as the most promising parent whereas “Granny Smith,” which has been phenotypically characterized as a regular phenotype (Lespinasse, 1977), did not cumulate the highest number of favorable alleles. “Coop17” was the most promising founder, estimated to be homozygous for the favorable allele at all QTLs. “Golden Delicious” and “Cox Orange,” widely used as founders in breeding programs (Noiton and Shelbourne, 1992), also carried a relatively high number of favorable alleles. Such reliable overview of the loci involved in bearing habit and estimation of the genotype at major loci is crucial for making relevant choices for breeding. The pedigree-based approach used herein takes the relationships between individuals into account by identity by descent and allows the transmission of alleles to be followed across a pedigree. This approach is particularly relevant for plant species in which varieties are tightly related to each other, which is the case in most crops (Soleimani et al., 2002). In addition, the relative importance of loci and the cumulative effects of small loci should not be overlooked. In this perspective, genomic selection models would be complementary to QTL analyses to evaluate the genetic value of individuals by summing allelic effects at each position of the genome (Kumar et al., 2013; Muranty et al., 2015).

## Author Contributions

EC planned and designed the research. AA and BG performed experiments and conducted fieldwork, JBD. performed statistical analyses, AA, Ev, and MB performed QTL detection. JBD, AA, BG, Ev, and EC interpreted the results and wrote the manuscript.

## Conflict of Interest Statement

Author Marco Bink is currently affiliated with Hendrix Genetics, however this work was completed whilst at a previous affiliation, Biometris, Wageningen University and Research centre, Droevendaalsesteeg. Therefore Marco Bink, and all other authors, have no competing interests to declare.

The reviewer PB and handling Editor declared their shared affiliation, and the handling Editor states that the process met the standards of a fair and objective review.

## Acknowledgments

The authors acknowledge fruitful discussions and advices from Y. Guédon and C. Trottier, S. Martinez, S. Hanteville, and Y. Holtz for their contribution in phenotyping and data analyses. The FruitBreedomics project N°. 265582 of the EU seventh Framework Programme (www.FruitBreedomics.com) is acknowledged for the usage of SNP, pedigree and genetic linkage map. The views expressed in this work are the sole responsibility of the authors and do not necessary reflect the views of the European Commission.

## Supplementary Material

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

## References

Aggelopoulou, A. D., Bochtis, D., Fountas, S., Swain, K. C., Gemtos, T. A., and Nanos, G. D. (2010). Yield prediction in apple orchards based on image processing. *Precision Agric.* 12, 448–456. doi: 10.1007/s11119-010-9187-0

Allard, A., Bink, M. C. A. M., Martinez, S., Kelner, J. J., Legave, J. M., Guardo, M., et al. (2016). Detecting QTLs and putative candidate genes involved in budbreak and flowering time in an apple multiparental population. *J. Exp. Bot.* 67, 2875–2888. doi: 10.1093/jxb/erw130

Antanaviciute, L., Fernández-Fernández, F., Jansen, J., Banchi, E., Evans, K. M., Viola, R., et al. (2012). Development of a dense SNP-based linkage map of an apple rootstock progeny using the *Malus* Infinium whole genome genotyping array. *BMC Genomics* 13:203. doi: 10.1186/1471-2164-13-203

Bates, D., Maechler, M., and Dai, B. (2011). *lme4: Linear Mixed-Effects Models Using S4 Classes. R Package Version 0.999375-28*. Available online at: http://lme4.r-forge.r-project.org (Accessed October 05, 2011).

Ben Sadok, I., Celton, J. M., Essalouh, L., El Aabidine, A. Z., Garcia, G., Martinez, S., et al. (2013). QTL mapping of flowering and fruiting traits in olive. *PLoS ONE* 8:e62831. doi: 10.1371/journal.pone.0062831

Bianco, L., Cestaro, A., Sargent, D. J., Banchi, E., Derdak, S., Di Guardo, M., et al. (2014). Development and validation of a 20K Single Nucleotide Polymorphism (SNP) whole genome genotyping array for apple (*Malus* × *domestica* Borkh). *PLoS ONE* 9:e110377. doi: 10.1371/journal.pone.0110377

Bink, M. C. A. M. (2002). On flexible finite polygenic models for multiple-trait evaluation. *Genet. Res.* 80, 245–256. doi: 10.1017/S0016672302005906

Bink, M. C. A. M., Boer, M. P., Braak, C. J. F., Jansen, J., Voorips, R. E., and van de Weg, E. (2008). Bayesian analysis of complex traits in pedigreed plant population. *Euphytica* 161, 85–96. doi: 10.1007/s10681-007-9516-1

Bink, M. C. A. M., Jansen, J., Madduri, M., Voorrips, R. E., Durel, C. E., Kouassi, A. B., et al. (2014). Bayesian QTL analyses using pedigreed families of an outcrossing species, with application to fruit firmness in apple. *Theor. Appl. Genet.* 127, 1073–1090. doi: 10.1007/s00122-014-2281-3

Bink, M. C. A. M., Uimari, P., Sillanpää, M. J., Janss, L. L. G., and Jansen, R. C. (2002). Multiple QTL mapping in related plant populations via pedigree-analysis approach. *Theor. Appl. Genet.* 104, 751–762. doi: 10.1007/s00122-001-0796-x

Breen, K. C., Tustin, D. S., Palmer, J. W., Hedderley, D. I., and Close, D. C. (2016). Effects of environment and floral intensity on fruit set behaviour and annual flowering in apple. *Sci. Hortic.* 210, 258–267. doi: 10.1016/j.scienta.2016.07.025

Burnham, K. P., and Anderson, D. R. (2002). *Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach, 2nd Edn*. New-York, NY: Springer-Verlag.

Bushakra, J. M., Sargent, D. J., Cabrera, A., Crowhurst, R., Lopez Girona, E., Velasco, R., et al. (2012). Rosaceae conserved orthologous set (RosCOS) markers as a tool to assess genome synteny between Malus and Fragaria. *Tree Genet. Genomes* 8, 643–658. doi: 10.1007/s11295-011-0450-y

Celton, J. M., Kelner, J. J., Martinez, S., Bechti, A., Khelifi Touhami, A., James, M. J., et al. (2014). Fruit self-thinning: a trait to consider for genetic improvement of apple tree. *PLoS ONE* 9:e91016. doi: 10.1371/journal.pone.0091016

Celton, J. M., Martinez, S., Jammes, M. J., Bechti, A., Salvi, S., Legave, J. M., et al. (2011). Deciphering the genetic determinism of bud phenology in apple progenies: a new insight into chilling and heat requirement effects on flowering dates and positional candidate genes. *New Phytol.* 192, 378–392. doi: 10.1111/j.1469-8137.2011.03823.x

Chagné, D., Crowhurst, R. N., Troggio, M., Davey, M. W., Gilmore, B., Lawley, C., et al. (2012). Genome-wide SNP detection, validation, and development of an 8K SNP array for apple. *PLoS ONE* 7:e31745. doi: 10.1371/journal.pone.0031745

Chan, B., and Cain, J. (1967). The effect of seed formation on subsequent flowering in apple. *J. Am. Soc. Hortic. Sci.* 91, 63–67.

Costes, E., and Guédon, Y. (2012). Deciphering the ontogeny of a sympodial tree. *Trees* 26, 865–879. doi: 10.1007/s00468-011-0661-8

Couranjou, J. (1983). Facteurs variétaux de l'alternance des arbres fruitiers. Résultats de quinze années de recherches chez le prunier domestique. *Fruits* 38, 705–728.

Dennis, F. G., and Neilsen, J. C. (1999). Physiological factors affecting biennial bearing in tree fruit: the role of seeds in apple. *HortTechnology* 9, 317–322.

Di Guardo, M., Micheletti, D., Bianco, L., Koehorst-van Putten, H. J., Longhi, S., Costa, F., et al. (2015). ASSIsT: an an automatic SNP ScorIng scorIng tool for in- and outbreeding species. *Bioinformatics* 31, 3873–3874. doi: 10.1093/bioinformatics/btv446

Di Pierro, E. A., Gianfranceschi, L., Di Guardo, M., Koehorst-van Putten, H. J., Kruisselbrink, J. W., Longhi, S., et al. (2016). A high-density, multi-parental SNP genetic map on apple validates a new mapping approach for outcrossing species. *Hortic. Res.* 3:16057. doi: 10.1038/hortres.2016.57

Durand, J. B., Guitton, B., Peyhardi, J., Holtz, Y., Guédon, Y., Trottier, C., et al. (2013). New insights for estimating the genetic value of segregating apple progenies for irregular bearing during the first years of tree production. *J. Exp. Bot.* 64, 5099–5113. doi: 10.1093/jxb/ert297

Dyk, M. M., van Soeker, M. K., Labuschagne, I. F., and Rees, D. J. G. (2010). Identification of a major QTL for time of initial vegetative budbreak in apple (*Malus* × *domestica* Borkh.). *Tree Genet. Genomes* 6, 489–502. doi: 10.1007/s11295-009-0266-1

Girona, J., Behboudian, M. H., Mata, M., Del Campo, J., and Marsal, J. (2010). Exploring six reduced irrigation options under water shortage for ‘Golden Smoothee’ apple: responses of yield components over three years. *Agric. Water Manage.* 98, 370–375. doi: 10.1016/j.agwat.2010.09.011

Gongal, A., Silwal, A., Amatya, S., Karkee, M., Zhang, Q., and Lewis, K. (2016). Apple crop-load estimation with over-the-row machine vision system. *Comput. Electron. Agric.* 120, 26–35. doi: 10.1016/j.compag.2015.10.022

Guitton, B., Kelner, J. J., Velasco, R., Gardiner, S. E., Chagné, D., and Costes, E. (2012). Genetic control of biennial bearing in apple. *J. Exp. Bot.* 63, 131–149. doi: 10.1093/jxb/err261

Hanke, M. V., Flachowsky, H., Peil, A., and Hättasch, C. (2007). No flower no fruit—genetic potentials to trigger flowering in fruit trees. *Genes Genomes Genomics* 1, 1–20.

Huff, A. (2001). A significance test for biennial bearing using data resampling. *J. Hortic. Sci. Biotechnol.* 76, 534–535.

Josse, J., and Husson, F. (2012). Selecting the number of components in principal component analysis using cross-validation approximations. *Comput. Stat. Data Anal.* 56, 1869–1879. doi: 10.1016/j.csda.2011.11.012

Kass, R. E., and Raftery, A. E. (1995). Bayes factors. *J. Am. Stat. Assoc.* 90, 773–795. doi: 10.1080/01621459.1995.10476572

Kumar, S., Garrick, D. J., Bink, M. C., Whitworth, C., Chagné, D., and Volz, R. K. (2013). Novel genomic approaches unravel genetic architecture of complex traits in apple. *BMC Genomics* 14:393. doi: 10.1186/1471-2164-14-393

Laurens, F., Aranzana, M. J., Arús, P., Bassi, D., Bonany, J., Corelli, L., et al. (2012). Review of Fruit genetics and breeding programmes and a new european initiative to increase fruit breeding efficiency. *Acta Hortic.* 929, 95–102. doi: 10.17660/ActaHortic.2012.929.12

Lauri, P. E., Combe, F., and Brun, L. (2014). Regular bearing in the apple – Architectural basis for an early diagnosis on the young tree. *Sci. Hortic.* 174, 10–16. doi: 10.1016/j.scienta.2014.05.001

Lauri, P. E., Crété, X., and Ferré, G. (2007). Centrifugal training in apple-appraisal of a two-year experiment on ‘Galaxy’ in Southeast France. *Acta Hort.* 732, 391–396. doi: 10.17660/ActaHortic.2007.732.59

Lauri, P. E., Gorza, O., Cochard, H., Martinez, S., Celton, J. M., Ripetti, V., et al. (2011). Genetic determinism of anatomical and hydraulic traits within an apple progeny. *Plant Cell Environ.* 34, 1276–1290. doi: 10.1111/j.1365-3040.2011.02328.x

Lauri, P. E., Terouanne, E., and Lespinasse, J. M. (1997). Relationship between the early development of apple fruiting branches and the regularity of bearing—An approach to the strategies of various cultivars. *J. Hortic. Sci.* 72, 519–530. doi: 10.1080/14620316.1997.11515539

Lespinasse, J. M. (1977). *La Conduite du Pommier (1ere Partie): Types de Fructification - Incidence sur la Conduite de l'arbre*. Paris: Invuflec.

Liu, W., Reif, J. C., Ranc, N., Porta, G. D., and Würschum, T. (2012). Comparison of biometrical approaches for QTL detection in multiple segregating families. *Theor. Appl. Genet.* 125, 987–998. doi: 10.1007/s00122-012-1889-4

Luan, T., Woolliams, J. A., Ødegård, J., Dolezal, M., Roman-Ponce, S. I., Bagnato, A., et al. (2012). The importance of identity-by-state information for the accuracy of genomic selection. *Genet. Sel. Evol.* 44, 1–7. doi: 10.1186/1297-9686-44-28

Molenberghs, G., and Verbeke, G. (2005). *Models for Discrete Longitudinal Data*. Springer Series in Statistics, New York, NY: Springer-Verlag.

Monselise, S. P., and Goldschmidt, E. E. (1982). “Alternate bearing in fruit trees,” in *Horticultural Reviews*, ed J. Janick (*Palgrave Macmillan*), 128–173.

Muranty, H., Troggio, M., Ben Sadok, I., Rifaï, M. A., Auwerkerken, A., Banchi, E., et al. (2015). Accuracy and responses of genomic selection on key traits in apple breeding. *Hortic. Res.* 2, 15060. doi: 10.1038/hortres.2015.60

Nagy, P. T., Szabo, Z., Nyeki, J., and Dussi, M. C. (2010). “Effects of frost-induced biennial bearing on nutrient availability and fruit disorders in an integrated apple orchard,” in *Xi International Symposium on Plant Bioregulators in Fruit Production*, ed G. Costa (Leuven: International Scociety of Horticultural Science), 753–757.

Noiton, D., and Shelbourne, C. J. A. (1992). Quantitative genetics in an apple breeding strategy. *Euphytica* 60, 213–219.

Pauly, L., Flajoulot, S., Garon, J., Julier, B., Béguier, V., and Barre, P. (2012). Detection of favorable alleles for plant height and crown rust tolerance in three connected populations of perennial ryegrass (*Lolium perenne* L.). *Theor. Appl. Genet.* 124, 1139–1153. doi: 10.1007/s00122-011-1775-5

Pearce, S. C., and Doberšek-Urbanc, S. (1967). The measurement of irregularity in growth and cropping. *J. Hortic. Sci.* 42, 295–305. doi: 10.1080/00221589.1967.11514216

Regnard, J. L., Segura, V., Merveille, N., Durel, C. E., and Costes, E. (2009). QTL analysis for QTL leaf gaz exchange in an apple progeny grown under atmospheric constraints. *Acta Hortic.* 814, 369–374. doi: 10.17660/ActaHortic.2009.814.60

Reddy, Y. T. N., Kurian, R. M., Ramachander, P. R., Singh, G., and Kohli, R. R. (2003). Long-term effects of rootstocks on growth and fruit yielding patterns of ‘Alphonso’ mango (*Mangifera indica* L.). *Sci. Hortic.* 97, 95–108. doi: 10.1016/S0304-4238(02)00025-0

Rosenstock, T. S., Rosa, U. A., Plant, R. E., and Brown, P. H. (2010). A reevaluation of alternate bearing in pistachio. *Sci. Hortic.* 124, 149–152. doi: 10.1016/j.scienta.2009.12.007

Samach, A., and Smith, H. M. (2013). Constraints to obtaining consistent annual yields in perennials. II: environment and fruit load affect induction of flowering. *Plant Sci.* 207, 168–176. doi: 10.1016/j.plantsci.2013.02.006

Segura, V., Cilas, C., Laurens, F., and Costes, E. (2006). Phenotyping progenies for complex architectural traits: a strategy for 1-year-old apple trees (*Malus* × *domestica* Borkh.). *Tree Genet. Genomes* 2, 140–151. doi: 10.1007/s11295-006-0037-1

Segura, V., Denancé, C., Durel, C.-E., and Costes, E. (2007). Wide range QTL analysis for complex architectural traits in a 1-year-old apple progeny. *Genome* 50, 159–171. doi: 10.1139/G07-002

Segura, V., Cilas, C., and Costes, E. (2008). Dissecting apple tree architecture into genetic, ontogenetic and environmental effects: mixed linear modelling of repeated spatial and temporal measures. *New Phytol*. 178, 302–314. doi: 10.1111/j.1469-8137.2007.02374.x

Segura, V., Durel, C.-E., and Costes, E. (2009). Dissecting apple tree architecture into genetic, ontogenetic and environmental effects: QTL mapping. *Tree Genet. Genomes* 5, 165–179. doi: 10.1007/s11295-008-0181-x

Smith, M. W., Shaw, R. G., Chapman, J. C., Owen-Turner, J., Slade Lee, L., Bruce McRae, K., et al. (2004). Long-term performance of ‘Ellendale’ mandarin on seven commercial rootstocks in sub-tropical Australia. *Sci. Hortic.* 102, 75–89. doi: 10.1016/j.scienta.2003.12.004

Soleimani, V. D., Baum, B. R., and Johnson, D. A. (2002). AFLP and pedigree-based genetic diversity estimates in modern cultivars of durum wheat [*Triticum turgidum* L. subsp. durum (Desf.) Husn.]. *Theor. Appl. Genet.* 104, 350–357. doi: 10.1007/s001220100714

Tustin, D. S., Dayatilake, G. A., Breen, K. C., and Oliver, M. J. (2012). Fruit set responses to changes in floral bud load-a new concept for crop load regulation. *Acta Hortic.* 932, 195–202. doi: 10.17660/ActaHortic.2012.932.28

van de Weg, E., Di Guardo, M., Koehorst-van Putten, H., Longhi, S., Noordijk, Y., Muranty, H., et al. (2013). *A Pipeline for Robust Marker Calling from Infinium SNP Arrays for Diploid Crops*. Available online at: https://hal.archives-ouvertes.fr/hajournalabbrevl-01209967

van de Weg, E., Voorrips, R. E., Finkers, R., Kodde, L. P., Meulenbroek, E. J., Jansen, J., et al. (2005). Pedigree genotyping: a new pedigree-based approach of QTL identification and allele mining by exploiting breeding material. *Acta Hortic.* 708, 483–488. doi: 10.17660/ActaHortic.2006.708.85

Velasco, R., Zharkikh, A., Affourtit, J., Dhingra, A., Cestaro, A., Kalyanaraman, A., et al. (2010). The genome of the domesticated apple (*Malus* × *domestica* Borkh.). *Nat. Genet.* 42, 833–839. doi: 10.1038/ng.654

Virlet, N., Costes, E., Martinez, S., Kelner, J. J., and Regnard, J. L. (2015). Multispectral airborne imagery in the field reveals genetic determinisms of morphological and transpiration traits of an apple tree hybrid population in response to water deficit. *J. Exp. Bot.* 18, 5453–5465. doi: 10.1093/jxb/erv355

Voorrips, R. E., Bink, M. C. A. M., Kruisselbrink, J. W., Koehorst-Van Putten, H. J. J., and van de Weg, W. E. (2016). PediHaplotyper: software for consistent assignment of SNP haplotypes in pedigrees. *Mol. Breed.* 36, 119. doi: 10.1007/s11032-016-0539-y

Keywords: bayes factor, biennial bearing, entropy, Malus × domestica, markov models, pedigree based analysis

Citation: Durand J-B, Allard A, Guitton B, van de Weg E, Bink MCAM and Costes E (2017) Predicting Flowering Behavior and Exploring Its Genetic Determinism in an Apple Multi-family Population Based on Statistical Indices and Simplified Phenotyping. *Front. Plant Sci*. 8:858. doi: 10.3389/fpls.2017.00858

Received: 15 December 2016; Accepted: 08 May 2017;

Published: 07 June 2017.

Edited by:

Katrin Kahlen, Hochschule Geisenheim University, GermanyReviewed by:

Ralf Uptmoor, University of Rostock, GermanyPeter Braun, Hochschule Geisenheim University, Germany

Copyright © 2017 Durand, Allard, Guitton, van de Weg, Bink and Costes. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Evelyne Costes, evelyne.costes@inra.fr

^{†}These authors have contributed equally to this work.