The present study was carried out on a population consisting of 265 DH lines derived from the F1 of a cross between homozygous parental lines P1 with good harvest-time reliability at higher temperatures and P2 that is less resistant to elevated temperatures but produces high-quality curds. More details about reliability in time to curd induction and harvest time of the PLs are provided by Hasan et al. (2016). A subset of 180 lines was used for model parameterization, while the validation sets included three subsets with 72 additional DH lines per se, 160 test hybrids from the parameterization set, and 34 test hybrids obtained from the 72 additional DH lines. Plant material was obtained from SYNGENTA Seeds B.V., Enkhuizen, Netherlands. The PLs were genotyped with a 20k B. oleracea Illumina Infinium iSelect single-nucleotide polymorphism (SNP) array (Matschegewski et al., 2015). A subset of SNPs of the array was selected based on polymorphisms between PLs and genome coverage to genotype the DH-population with 176 SNP markers. The genetic map was constructed using MapMaker 3 and the Haldane function. The total map length was 891.2 cM spanning over nine chromosomes. Average distance between SNP markers was 5.3 cM. Physical marker positions are available from Hasan et al. (2016).

Model parameterization is based on controlled temperature treatments described in Hasan et al. (2016). Six trials were carried out in greenhouses with mean temperatures of 11.8, 15.5, 17.3, 19.0, 21.4, and 27.0°C and one in climate chambers with a mean air temperature of 26.5°C. Seeds of all DH lines of the parameterization set were sown in seedling trays and raised in the greenhouse at 22°C until all plants had two to three visible leaves. Afterward, seedlings were planted into 3 l pots, transferred to the temperature treatments in greenhouses or climate chambers, and arranged in randomized complete block designs. Plants were grown at a photoperiod of 16 h, providing 200 μmol m^-2 s^-1 additional light by Phillips SON-T Agro lamps when day length was below 16 h or natural radiation during the day below 5 klux. Plants were fertigated daily with 0.5 g l^-1 Scotts Universal solutions. All trials were terminated after 120 days.

Time to visible curd induction was recorded daily. A curd diameter of 1 cm was used as threshold for considering visible curd induction to have been occurred. Number of leaves larger than 1 cm was counted twice a week. Leaf appearance rate (LAR) was computed as the slope of the regression of leaf numbers on growing degree-days assuming a base temperature of zero.

For model validation, a field trial was carried out with 72 lines in 2013 in Rostock, Germany. Twelve seeds per line were sown on May 6 into seedling trays and raised in the greenhouse. Eight plants of each line were planted into the field on June 5 when most plants had two to three leaves. Field trials were carried out in a randomized complete block design with two replications and four plants per replication. Seedlings were planted in rows with distances of 50 cm between and within rows. A border row was planted to provide similar conditions for all plants. Irrigation was applied as required and the plot was covered with a net to protect plants during the first weeks from cabbage fly (Delia radicum) and other pests. The trial ended on August 16. Time to curd induction was scored twice a week. Beginning of curd formation was defined as curd diameter ≥ 1 cm.

Additional validation trials with similar growing conditions were carried out in 2013 in Ruthe, Germany and in 2011 and 2012 in Zeewolde, Netherlands. In Ruthe, seeds of the validation set were sown on May 14 and eight plants per line were transplanted to the field on June 13 when they had two to three visible leaves. Zeewolde field trials were carried out on the whole DH-population and transplanting took place on June 14, 2011 and May 24, 2012. Mean air temperatures during the experiments were 17.9°C in Rostock, 19.0°C in Ruthe, and 17.4 and 17.2°C in Zeewolde 2011 and 2012, respectively. Only 44 DH lines were used in Zeewolde 2011, whilst all other trials comprised the whole DH per se validation set.

Further trials were conducted on the 160 F1 test-hybrids derived from the parameterization set and the 34 hybrids in Ócsa, Hungary, where plants were transplanted on July 7, 2011 and July 6, 2012 and in Zeewolde, where transplanting took place on June 13, 2012. Mean air temperatures in Ócsa were 21.4 and 22°C, respectively; mean air temperature in Zeewolde was 17.4°C. Harvest time was recorded in all Ósca and Zeewolde trials. Sowing was assumed to have taken place 30 days before transplanting and duration from visible curd induction to harvest was assumed to be 30 days (Table 3).

The cauliflower phenology model was divided into a juvenile and a vernalization phase. During the juvenile phase, plants were assumed to be insensitive to temperature, i.e., development followed the temperature-sum rule. Although it is known that genotypic variation in the number of initiated or visible leaves at juvenile-to-adult transition exists (Hand and Atherton, 1987), juvenility was assumed to end for all genotypes after 7 leaves were visible. In contrast, LAR and phyllochron (P, d), the reciprocal of LAR, was assumed to be genotype specific. The genotype specific temperature sum for juvenile-to-adult transition (TS_f) was calculated as follows:

The model assumed genetic variation in temperature response from end of juvenility to visible curd induction (curd diameter > 1 cm). Accelerating or retarding temperature effects on phenological development may end at floral transition, which, however, was not observed. Instead, the relation between thermal time to visible curd induction in seven different treatments and temperature was estimated by linear regression, where the intercept T₀ (°Cd) describes the extrapolated (theoretical) minimum thermal time to curd induction and the slope S (°Cd °C^-1) the temperature sensitivity, i.e., a genotype with S = 0 is insensitive while a genotype with a large positive S shows strong vernalization response. The daily development rate k_ij (d^-1) for the i^th day and j^th line was calculated based on daily mean air temperatures T_i as follows:

k_ij was added up until the cumulative k_j became ≥ 1, i.e., when the adult vegetative phase was completed and curd development started.

A linkage map was constructed with MapMaker 3 by simultaneous multipoint-analysis using the Haldane function (Lander et al., 1987; Lincoln et al., 1993). On basis of the results of parameterization trials, QTL analyses were conducted on LAR, S, and T₀ using PlabQTL 1.2 (Utz and Melchinger, 1996). Composite interval mapping (CIM) based on multiple regression with co-factors was applied (Haley and Knott, 1992; Jansen and Stam, 1994; Zeng, 1994; Miller, 2002). The LOD threshold to define a significant QTL was 2.5. The linkage map was drawn with MapChart 2.3 (Voorrips, 2002).

Additive QTL effects were assigned to all detected loci and QTL based parameters for T₀, S, and LAR were calculated using the following equation:

where y_j is the estimated value of parameter y of the j^th line, m is the population mean for the parameter, e_i is the additive effect of the i^th QTL, and g_ij is the allelic state of the i^th QTL in line j. g_ij can be -1 or +1.

In addition to linkage mapping, marker effects of the three parameters were estimated by genomic selection using the ridge-regression best linear unbiased prediction (rrBLUP) package for R 3.2.4. with the following restricted mixed model (REML) solution:

where y is the vector of phenotypic observations, X is the matrix of fixed effects β, Z is the matrix of random effects υ, and ε is the residual effects matrix (Endelman, 2011). rrBLUP was carried out on the three model parameters LAR, S, and T₀. Effects were assigned to all markers and the sum of marker effects was used to calculate additive genomic estimated breeding values (GEBV) of the individuals of the validation sets. GEBVs for LAR, S, and T₀ were implemented into the phenology model as parameter input values to simulate heading date of the individuals in field trials.

Thermal time from transplanting to curd induction under different constant temperatures was measured in the parameterization trials. In general, the higher the temperature, the higher was the temperature sum required for curd induction. Population means ranged from 753 to 1590°Cd (Table 1).

Mean LAR of individual genotypes averaged across all parameterization trials ranged from 0.016 to 0.029. The population mean was 0.022 with a standard deviation of ±0.002. Mean S was 50.0 (±20.5) and ranged from 4.2 to 106.0. Population mean of T₀ was 631°Cd ± 296°Cd (Table 2). Mean R² of the relation between temperature and thermal time to curd induction was 0.65 (±0.13), maximum R² was 0.87 and minimum was 0.05.

In the validation trials the thermal time to curd induction was measured for the parameterization set, the validation set and F1 hybrids of both sets. Mean thermal time to curd induction ranged from 1276.3 to 1597.7°Cd across all trials and the different subsets (Table 3).

QTL analyses were carried out on model parameters LAR, S, and T₀. In total, ten QTL were detected for LAR. Three of them were located on chromosome O1 and two on O5 (Figure 1). Five QTL were identified for S. Both S QTL on O6 co-localized with T₀ QTL, while one of them was also found to co-localize with an LAR QTL. The third QTL for S on O9 was located close to a LAR QTL. Genetic positions, LOD scores, and additive effects of all QTL are summarized in Table 4. LAR QTL on O9 and O5 revealed highest LOD scores of 7.04 and 6.62, respectively. LOD scores of the two QTL for S on O6 were found to be highest with 8.0 and 4.6.

Genome-wide marker effects for the model parameters LAR, S, and T₀ were estimated by rrBLUP and are displayed in Figure 2. Large marker effects on O6 support linkage-mapping results, which indicated a promising hotspot region on O6.

To get an idea about the power of the phenological model, curd induction of the parameterization set was simulated with both QTL-effects and original parameters. A very high prediction accuracy of R² = 0.98 was observed applying original parameters to predict time to curd induction of the parameterization trials. Figures 3A,C show means across seven temperature treatments. R² decreased to 0.44 if original parameters were replaced by QTL-based estimates. A similar situation with an overall reduced R² was observed if mean time to curd induction of field data from Zeewolde 2011 and 2012 was simulated (Figures 3B,C). R² diminished to 0.47 if mean curd induction of F1 hybrids derived from the parameterization set in field trials in Ócsa 2011, 2012 and Zeewolde 2012 was modeled using original parameters (data not shown).

For model evaluation, QTL effects were used as parameter inputs to predict time to curd induction of the validation set. Figure 4 shows predicted vs. observed time to curd induction of different field trials. In Rostock, mean observed time to curd induction was 83.5 days (±11.5), while predicted time to curd induction was 79.1 days (±4.9). The coefficient of determination was R² = 0.45 (Figure 4A). With R² = 0.47, prediction accuracy was slightly higher in Ruthe (Figure 4B), while the model performed best in the Zeewolde trial in 2011 (R² = 0.51, Figure 4C). Interestingly, prediction of independent lines resulted in higher R² values if compared to the parameterization set. However, the slope was 1.08 for the parameterization set (Figure 3D) but ranged between 1.22 and 1.54 when curding time of independent lines was modeled (Figure 4).

As an alternative method, genome-wide marker effects were estimated. The GS model led to higher prediction accuracies in all field trials if applied on the validation set (Figure 5). R² increased on average by 10.3%. However, a relatively strong bias was observed in the Zeewolde trials (Figures 5C,D).

Applying the QTL-based model on test hybrids developed from the parameterization set resulted in moderate prediction accuracies (R² = 0.34, Figure 6A). If QTL effects were replaced by genome-wide marker effects, R² increased to 0.50 (Figure 6C). The QTL based simulation of time to curd induction of F1-hybrids derived from the validation set revealed also moderately good results (R² = 0.41, Figure 6B), which were outperformed by the GS based model (R² = 0.48; Figure 6D).

QTL were found for all examined parameters. QTL for S and T₀ were almost congruent in their positions, which was not surprising, as S and T₀ were highly correlated (r = -0.95) and lines with a steep slope do have a low intercept and vice versa. Similar to the present study, Uptmoor et al. (2008) found QTL for ‘slope’ on chromosomes O9 and O6 in a B. oleracea var. alboglabra × B. oleracea var. italica cross. QTL for ‘slope’ were not detected on chromosomes O4 or O1 but an additional QTL was identified on O3, which was not confirmed in the present study. Uptmoor et al. (2008) detected four QTL for the parameter ‘intercept.’ One was also located O6 but three further QTL on O3, O4, and O5 were not verified in this study. Even if positions are not completely comparable since both studies used different markers and linkage maps, clear parallels between the detected QTL are evident.

At the positions of three of the five QTL for slope, across environment QTL with significant QTL × environment interactions were detected in an earlier study using the same population and phenotype data (Hasan et al., 2016), indicating that input parameters of ecophysiological models are suitable traits for QTL detection. In contrast to the statistical interaction of a QTL with its environments (i.e., different effect sizes in different environments), model parameters display the trait’s response to an environmental factor allowing precise estimations of a new genotype in untested environments (Onogi et al., 2016). The QTL regions on O6 harbor BoAP1-c and BoAP1-a, the latter of which was suggested to interact with temperature (Labate et al., 2006). The QTL region on O9 harbors an FLC paralog and a FRIGIDA LIKE (FRL) ortholog (Hasan et al., 2016).

Usually, a relatively large sample size is required for QTL analyses and several hundreds of lines are recommended for good results (Hyne et al., 1995; Beavis, 1998; Lynch and Walsh, 1998). In general, the higher the number of genotypes used for analyses, the more accurate are the results. QTL with relatively small effects, so called minor QTL, are detected only in large populations, while major QTL can still be found with fewer genotypes (Vales et al., 2005; Slate, 2005), which, however, leads to an underestimation of the total number of QTL, while QTL effects are overestimated since a part of the variation is due to undiscovered minor QTL (Beavis, 1994; Xu, 2003). All significant QTL were tested for QTLxQTL interactions by pairwise linear regression analysis. The only significant interaction was detected for the parameter slope (QTL at 114 cM on O4 and at 12 cM on O6). However, including the interaction effect into the model improved estimations for the slope only slightly and had no positive impact on the simulation studies (data not shown). Hill et al. (2008) found that additive variance accounts often for close to 100% of the total genetic variance, which might be the reason for the missing positive impact of QTLxQTL interactions on the QTL-based simulation model of the present study.

Genomic selection is a useful tool to overcome these limitations as effects are assigned to all markers. GS was used as an alternative method for model parameterization and regions with large marker effects on O4 and O6 reflect findings from linkage mapping quite well; but in contrast to the QTL model, GS still accounts for genome regions with minor effects.

The original model was able to predict time to curd induction by means of chosen parameters in a satisfactory manner. The possibility to simulate cauliflower development from temperature data only has already been shown (Olesen and Grevsen, 2000). The QTL-based approach was rather successful as well, although predictions based on original parameters were more precise. Assumptions made by Reymond et al. (2003) and Uptmoor et al. (2008) that replacement of model parameters by QTL effects would lead to higher prediction accuracies due to the reduction of random errors were not confirmed in the present study. While the prediction of hybrids derived from the DH population was possible, prediction accuracies were not as high as those of DH lines per se.

Considering simulation studies carried out on independent genotypes of the cross, results adequately reflected expectations if compared to similar studies without independent genotypes (Quilot et al., 2005; Yin et al., 2005a; White et al., 2008; Uptmoor et al., 2012). However, Bogard et al. (2014) attained relatively high prediction accuracies (mean R² = 0.58) for predicting heading date of independent bread wheat genotypes using a common phenology model and genome-wide association studies for QTL detection. Less precise predictions of the present study may result from existing variation in juvenile-to-adult transition, for which the curd-induction model did not account.

Unexpectedly, the genome-based prediction of independent DH lines exceeded the accuracy of the simulations of lines used for parameterization and even the prediction of F1 hybrids resulted in a slightly higher R². If subsamples of a population are large enough, different random subsamples will provide similar results regarding QTL positions and effects, hence, prediction accuracies for independent lines are likely to be on the same level as for lines used as parameterization set. Therefore, higher accuracies for independent lines in this study can only be explained by chance. In fact only R² values were higher, while slopes > 1 illustrate that the variability in curding time within the validation set was overestimated.

Simulations using genome-wide marker effects led to more precise predictions in comparison to QTL-based simulations, which held true for both independent DH lines of the validation set and test hybrids. The major reason might be explained by missing effects of minor loci in the QTL model. Comparing different models to predict hybrid performance in maize, Guo et al. (2013) found the genome-wide prediction approach being more accurate than predictions based on QTL effects.

In the present study, a simple approach with only three model parameters was used to predict time to curd induction of cauliflower lines and hybrids. Combining GS and ecophysiological models could reduce efforts spent in field trials to test germplasm across wide ranges of different environments in silico since GS models are able to predict phenotypes only in relatively similar environments while integrated models account for both genotypic and environmental variations in equal measure. Even though model parameterization requires large phenotyping efforts under controlled conditions, more precise predictions may be attained from more complex models.

The uncertainty about juvenile-to-adult transition might be considered as main error source impairing prediction accuracies. Juvenility is assumed to be independent from temperature and to end after a certain number of leaves are initiated (Salter and James, 1974; Hand and Atherton, 1987). Since large variations in juvenile-to-adult transition are present in cauliflower cultivars (Booij and Struik, 1990), certain variation may exist in the examined population as well. Genotype specific estimations of leaf numbers at the end of juvenility should be made in further experiments via reciprocal transfers to improve the model’s accuracy by integrating the onset of sensitivity to temperature.

In the present model, LAR was assumed to be constant throughout the vegetative phase. However, it is known that LAR changes with time and more sophisticated models already took this into account (Kage and Stützel, 1999). Similar approaches are conceivable to further improve simulation accuracies of integrated modeling approaches.

While the model of the present study simulates the development of cauliflower until visible curd induction, other phenology models already implemented curd development to predict harvest time (Pearson et al., 1994; Olesen and Grevsen, 2000). The practical use of a model predicting harvest time instead of time of curd induction is undoubtedly higher since accurate harvest-time predictabilities in diverse environments are a great challenge not only for growers but also for breeders.