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).

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.

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).

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 Ent¯g(entropy based on frequencies) and Ent¯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 μ₀₀ = 0), φ_t is the fixed effect of year t (with reference φ₂₀₀₆ = 0) treated as a qualitative variable, variables θ_{g, m} are independent random interactions between genotype g and memory m with common variance τθ2, variables η_{g, t} are independent random interactions between genotype g and year t with common variance τη2, and variables ζ_{g, r} are independent replication-specific random effects with common variance τζ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.

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, Ent¯gand Ent¯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.

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, Ent¯gand Ent¯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 Ent¯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.

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, Ent¯g, Ent¯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,σa2I) and e~N(0,σe2I). σa2 and σe2 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(μ,a,λ,σa2,σe2|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).

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 τθ2=τη2=τζ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.

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 μ₁₀ was higher than μ₀₁ and μ₁₁), 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.

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 Ent¯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 Ent¯g, and 0.64 using PCs. Using Ent¯glmm,g instead of (or in addition to) Ent¯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 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 Ent¯g, and 35% with PCs. Genotypes with unknown values of BBI_res_norm_ax, γ^ax or Ent¯gwere excluded, thus keeping 115 genotypes in SG. As previously, the best prediction was achieved with three PCs. Both Ent¯glmm,g and Ent¯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.

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}, Ent¯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 Ent¯g only)

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 Ent¯g were moderately correlated with BBI_res_norm (negatively) and with γ^ax (positively), while Ent¯glmm,g was poorly correlated and θ_{g, 10} uncorrelated (at level 0.05) with the tree-scale indices.

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.

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 Ent¯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 Ent¯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 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).

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 Ent¯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).

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 Ent¯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).

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 Ent¯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.

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 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.