Introduction

The study of developmental processes in plants requires the evaluation of traits over time taking into account the continuous nature of development. Quantitative traits normally involved in these processes require assessments at several time points during the life cycle. However, conventional experiments include evaluations at a single fixed time point during the growing season. Quantitative trait loci (QTL) analysis using the data collected from those experiments give just an impression of the loci affecting the trait at a particular developmental stage. The understanding of the genetic basis controlling quantitative traits improves when the evolution of the trait during the life cycle and the developmental pattern are considered. In addition, environmental factors affecting crop development, such as temperature and photoperiod, must be also taken into account. Therefore, growth models flexible enough to be adapted to the usually non-linear trait responses over time have to be implemented. Examples of commonly used models in biology are the exponential growth models and the family of s-shaped curves (Schnute 1981). These models have the advantage of describing the development of the trait in terms of curve parameters with a biological interpretation. The differences in growth trajectories between genotypes are reflected by genotype-specific curve parameters. Growth models and QTL analysis can be combined by modeling growth curve parameters in terms of QTLs in one-step (Ma et al. 2002; Malosetti et al. 2006) or two-step model approaches (Reymond et al. 2003; Yin et al. 1995). In the former, logistic growth curves and QTL mapping are combined, modeling growth parameters as a function of the QTL genotypes described by molecular marker scores. In the latter, first genotype-specific parameters are estimated from the growth curves, thereafter they are used as phenotypic traits in a conventional QTL analysis. In both approaches, the accurate estimation of the curve parameters plays an important role to find meaningful results.

In our study, haulm senescence was assessed at several time points during the growing season on a diploid potato mapping population evaluated for two consecutive years in Finland (2004 and 2005). We are presenting here an alternative two-step approach: modeling first in a flexible way the curve trajectories with a smoothed generalized linear model (GLM) procedure. We used the collected time series data to model the development of senescence in terms of time and temperature. To allow a better comparability of the trait across the years, we converted the calendar days after planting into beta-thermal days after planting (TAP) using a non-linear temperature effect function described by Khan et al. (2011). Once the senescence trajectory was defined, the characteristics of the curve-such as onset, maximum and average progression rate and inflection point- were estimated from the first and second derivative of the curve. These characteristics were used as phenotypic traits describing senescence development in the QTL analysis. In this way, the characterization of the smoothed GLM curves allowed us to map the genetic basis of the senescence process giving a meaningful biological interpretation to the results.

Materials and methods

Description of the C×E potato population

The evaluation of haulm senescence was done on a diploid backcross potato population (C×E) composed of 250 genotypes. The C×E population is the result of a cross between two diploid potato clones. The female parent, the C clone US-W5337.3 (Hanneman and Peloquin 1967) is a hybrid between Solanum phureja (PI225696) and S. tuberosum dihaploid US-W42. The male parent, the E clone (Jacobsen 1980), is a hybrid between VH34211 (a S. vernei–S. tuberosum backcross) and the C clone itself. S. phureja is characterized by a lack of dormancy, early maturity and short day tuberization induction. In contrast, S. tuberosum varieties have long dormancy, variable maturity and long day tuberization induction (Hawkes 1990; Ewing and Struik 1992). The different day length requirements for tuberization of S. phurja and S. tuberosum has provided a source of genetic variation in the C×E population making it suitable for the study of developmental processes under different photoperiods.

Experimental design and description of the haulm senescence evaluation

The very long day conditions of northern Finland were used to conduct an experiment in two consecutive years (2004 and 2005) in the experimental field of AgriFood Research Finland (MTT) at the North Ostrobothnia Research Station in Ruukki (64°42′N, 25°00′E). The experiment performed in 2004 was previously described by Zaban et al. (2006). Sets of 197 and 222 genotypes of the C×E population were planted in plots of three plants per genotype, where the genotypes were randomized, on June 1st 2004 and May 16th 2005, respectively. We acknowledge that in these trials there was no real replication of genotypes, but only pseudo-replications. Each plant was evaluated at several time points spaced in intervals of 3–7 days. In 2004 the observation period was 52 days with evaluations at 77, 81, 86, 91, 95, 100, 105, 109, 116, 123, 129 days after planting (DAP). In the second year the observation period was shorter (31 days) and haulm senescence was evaluated at 84, 88, 93, 98, 101, 106, 109, 112 and 115 DAP. The process of senescence was defined as the period between the last observation at which the plant was entirely green and the first date at which the plant was dead (Celis-Gamboa 2002). The progress of the trait was measured using the scale described by Celis-Gamboa et al. (2003) in which 1 = green plant; 2 = upper leaves with the first signs of yellowing (light green); 3 = yellow leaves; 4 = 25% of haulm tissue brown; 5 = 50% of haulm tissue brown; 6 = more than 75% of haulm tissue brown; 7 = dead plant.

Description of Beta-thermal time estimation

Crop development is mainly affected by temperature and can be modified by other factors such as photoperiod (Hodges 1991). Previous studies have shown that potato growth is influenced by temperature, where warmer temperatures favor vegetative growth accelerating reproductive and vegetative development (Haun 1975; Benoit et al. 1986; Struik and Ewing 1995). Whereas, lower temperatures facilitate tuber growth (Marinus and Boadlaender 1975). The effect of temperature on crop development rate is often described by using a thermal time concept. Various non-linear models have been developed to describe temperature response of developmental processes in plants (Johnson and Thornley 1985; Gao et al. 1992; Yin et al. 1995). In our study, the daily contribution of temperature to plant development in 2004 was different from 2005 due to the fluctuations in daily air temperature. A summary of the meteorological conditions reported by the station in Ruukki (Finland) during the growing season in the two experimental years is presented in Table 1. The growing season defined as the period between the last killing frost of spring and the first killing frost of autumn (Allaby 1998) was defined in 2004 from April 16 to November 10 and in 2005 from May 3 to November 14. To allow for a better comparability of senescence across the years, we converted the calendar DAP into beta-thermal days after planting (TAP) using the non-linear temperature effect function g(T) described by Khan et al. (2011):

$$ g(T) = \left[ {\left( {\frac{{T_{c} - T}}{{T_{c} - T_{o} }}} \right)\left( {\frac{{T - T_{b} }}{{T_{o} - T_{b} }}} \right)^{{\frac{{T_{o} - T_{b} }}{{T_{c} - T_{o} }}}} } \right]^{{c_{t} }} $$
(1)

where the three main temperatures, base, optimum and ceiling temperature for phenological development of potato, were defined as T b  = 5.5°C, T o  = 23.4°C and T c  = 34.6°C, respectively. The temperature response curvature coefficient was estimated as c t  = 1.7 according to Khan et al. (2011). Because the function g(T) is non-linear and temperature fluctuates daily, g(T) was estimated using the average daily air temperature to obtain the daily value. The cumulative TAP, combining temperature and time (beta-thermal time, BTT), was the scale of the x-axis used to compare senescence development in the two experimental years.

Table 1 Summary of meteorological conditions reported in Ruukki, Finland during the growing season in 2004 and 2005

Description of the smoothed generalized linear model procedure

To estimate the senescence curve non-parametrically, a variant of penalized B-splines, P-splines (Eilers and Marx 1996), was used. Then we interpreted a senescence value y as having observed y “successes” in n binomial trials and we estimated a smooth curve for the probability p of a “success”. To guarantee that 0 < p < 1, we worked with the linear predictor. In this case it is the logit, η = log(p/(1−p)), which has no restrictions on its range. This is standard practice in logistic regression and generalized linear models (Nelder and Wedderburn 1972). Smooth logistic regression with P-splines was explained by Eilers and Marx (1996). They also showed that automatic interpolation was obtained. We exploited this property by working with zero degree B-splines including a knot on each beta-thermal day (TAP units) of the domain we study. The penalized log-likelihood in this situation is

$$ l^{*} = \sum\limits_{i} w_{i} [y_{i} \log p_{i} + (n_{i} - y_{i} ) \log ({1- p_{i}})+ \log \left( \begin {array}{ll}{n_{i}}\\ {y_{i}} \end{array} \right) ]-\frac {\lambda} {2}\sum\limits_{i} {(\Delta ^{2}\eta _{i} )^{2}} - {\frac {\kappa} {2}}{\sum\limits_{i} {\eta_{i}^{2}}} $$
(2)

with \( p_{i} = 1/\left( {1 + e^{ - \eta } i} \right). \)The weight is determined by \( w_{i} = np_{i} \left( {1 - p_{i} } \right) \) following GLM methodology for our situation. The second term in the equation above is the penalty on the linear predictor η. The parameter λ tunes its weight: the larger λ, the smoother the result. We used the operator notation to indicate repeated differences: \( \Updelta^{2} \eta_{i} = \eta_{i} - 2\eta_{i - 1}\,+\,\eta_{i - 2} = \left( {\eta_{i} - \eta_{i - 1} } \right) - \left( {\eta_{i - 1} - \eta_{i - 2} } \right).\)

One can show that for very large λ the estimated curve for p approaches the logistic curve. Furthermore we used a ridge penalty on η tuned by the parameter κ to avoid numerical instabilities, when η becomes too large. To estimate η, the iterative weighted linear regression algorithm (Nelder and Wedderburn 1972) was modified, to account for the penalties.

The fitting procedure resulted in a smooth curve for the development of haulm senescence over time for each genotype. This flexible functional form using splines allowed to model different shapes of these curves (e.g., for early and late genotypes) within the same framework.

Characteristics estimated from the senescence curves

Once the senescence curves were fitted, some characteristics of the curves describing the ageing process were estimated to allow the study of senescence as a continuous trait changing in time. The characteristics estimated from the curve have a meaningful biological interpretation for the senescence development and facilitate the understanding of the QTL mapping results. We calculated the first differences on the curve values (first derivative) as a proxy of the slope of the curve. The mean slope or mean progression rate (mprate) is reflecting the average rate of change of the senescence curve during the whole observation period giving an idea of how fast a genotype is experiencing the senescence process. The maximum slope (prate) explains the maximum rate of change of the senescence curve when the plant has completed half of the senescence process (between 4 and 5 in the senescence scale). The higher the value, the faster full senescence is reached during the observation period. The inflection point (ipoint) reflects the point in time when half of the senescence process has been reached and the trajectory curve change from convex to concave shape indicating the beginning of the final stage. Additional characteristics were deduced from the second differences of the curve values (second derivative), including the maximum and minimum change of the slope, which are interpreted as the onset and the end of the senescence process. The onset is indicating the beginning of the senescence process in terms of time or more accurately in our case in terms of BTT. The lower the value, the earlier the senescence process starts.

Repeatability estimation

Data from the three plants within a plot for individual genotypes were used for a repeatability estimation, calculated for each year. The repeatability was defined by the ratio of genetic (genotypic) to phenotypic variance. Phenotypic variance was equal to the sum of the genetic variance and a third times the between-plants-within-plot variance.

G×E interactions

The estimation of genotype by year interactions included 186 genotypes evaluated in the two consecutive years using a two-way analysis of variance with genotype and environment fixed in the model. The curve characteristics (mprate, prate, ipoint and onset) were used as response variables, y, and the genotype (G), year (E) and Genotype × Year (G × E) interaction as the explanatory variables in the model, which also includes an error term (based on the variation between the three plants per genotype):

$$ y = {\it{G}} + {\it{E}} + {\it{G}} \times {\it{E}} + {\it{error}}.$$
(3)

Description of the genetic map and molecular data

The C×E population was genotyped using amplified fragment length polymorphism, AFLP (Celis-Gamboa 2002), simple sequence repeat, SSR, and cleavage amplified polymorphism, CAPS (Werij et al. 2007). Our study included a subset of dominant and co-dominant markers with the expected segregation ratios 1:1 and 1:1:1:1, respectively. JoinMap 4 (Van Ooijen 2006) was used to construct the C and E parental maps using 164 and 198 markers, respectively, each with 12 linkage groups (LG) previously reported by Celis-Gamboa (2002). The C map consisted of 135 markers spanning 917.4 cM. Four of the 12 LG were split in two sub-groups and one was split in 3 sub-groups due to the large distances between adjacent markers (more than 30 cM). The E map consisted of 132 markers spanning 629.8 cM and 2 of the 12 LG were split in 2 sub-groups. Since the maternal (C) and paternal (E) maps were not integrated due to the expected differences in the recombination frequencies between the two parents (with genetic background from two different Solanum species), the separated parental maps were used to perform the QTL analysis. Co-dominant markers present in both parental maps (C and E) were used to identify the same LG in the two maps. The assignment of linkage groups was done using as a reference the maps of Celis-Gamboa (2002), each LG is preceded by the letter C or E according to the parental map, followed by the LG number.

QTL analysis using curve characteristics and single time points

The detection of QTLs, estimation of QTL main effects and their interactions was done in two steps. In the first step, the detection of QTLs was performed separately for each type of phenotypic data available. In the analysis per time point the senescence scores, from 1 to 7, were used in the QTL mapping. In the analysis of the characteristics estimated from the curves (mprate, prate, ipoint and onset), each characteristic was used as a quantitative trait. The use of model parameters in QTL mapping has been introduced in functional mapping for the analysis of time-series traits. It incorporates biological principles (defined with mathematical functions) into the framework for QTL mapping (Wu et al. 2003; Wu and Lin 2006).

A nonparametric QTL analysis was performed for the time point data using the rank sum test of Kruskal–Wallis (KW). The curve characteristics were analyzed using both KW and interval mapping (IM). Both procedures are available in MapQTL 6 (Van Ooijen 2009). For the curve characteristics, QTLs detected by KW and IM were the same. For transparent comparison of QTL analysis on time points and curve characteristics, in this paper we present only results from KW.

QTL mapping was performed separately on the parental maps and the criterion for detecting QTLs was set at a significance level of P ≤ 0.005 (Van Ooijen 2006). To identify the detected QTLs, the origin of the map (C and E parent) and the linkage group are indicated in the QTL name using italic letters. The QTLs were named according to the parental map and the linkage group in which they were identified. As an example, the most significant QTL on linkage group 4 in the C parent was called C4 and in the case of common markers detected in both parents, like E5 and C5, the QTL was called Ch5 (corresponding to a QTL on chromosome 5).

QTL × QTL interactions

IM was used for detecting main effect QTLs for the curve characteristics mprate, prate, ipoint and onset. The presence of epistatic interactions was investigated per year by fitting linear models (Genstat release 13.2, Payne et al. 2010) containing pairwise interactions between those QTLs for whom earlier main effects were detected for any of the four curve characteristics. A full model was defined with each curve characteristic as response variable, y, and the main effect QTLs together with all possible QTL × QTL interactions as fixed explanatory variables.

$$ y = {QTLs} + {QTL} \times {QTL \; interactions} + {error} $$
(4)

From the full model (4), for each curve characteristic a subset of model terms was created. This list contained the significant main effects QTLs, all significant epistatic interactions, plus non-significant main effect QTLs whenever these QTLs were involved in significant epistatic interactions. A final model was fitted for each curve characteristic including only the earlier selected terms.

Results

Curve fitting

The transformation of normal calendar days after planting into TAP and their cumulative values generated a new scale (BTT) to measure the senescence development in terms of time and accumulated temperature. BTT went from 33.3 to 41.85 and from 37.1 to 50.8 during the observation period in 2004 (77–129 DAP) and 2005 (84–115 DAP), respectively. This scale was used as the x-axis to fit the curves and visualize the progression of haulm senescence scored from 1 to 7. The use of the smoothed GLM procedure in our study allowed flexible modeling of different curve shapes. Different types of curves were observed in the C×E population according to the maturity type of each genotype (Fig. 1a). The description of the five categories (very early, early, intermediate, late and very late) according to the duration of the plant cycle was reported by Celis-Gamboa (2002). In early genotypes we observed a logistic s-shape in the curve trajectory while in the late genotypes we had only the exponential part of the curves often corresponding to completion of only about half of the senescence process.

Fig. 1
figure 1

Characterization of senescence development in the C×E population represented by three genotypes with different maturity type: CE78 (early: red-top line), CE47 (intermediate: blue-middle line) and CE697 (late: yellow-lowest line). a Fitted curves showing the senescence development. b First derivative of the curves in which the maximum slope represents the progression rate and the dashed lines represent the time when half of the senescence process has been completed (inflection point) in each genotype. c Second derivative of the curves in which the maximum second difference is interpreted as the onset of senescence and it is indicated with a dashed line in each genotype. BTT beta-thermal time representing cumulative TAP beta-thermal days after planting

After fitting the curves for each genotype, several characteristics describing the ageing process were estimated to be used in the QTL analysis. The first differences on the curves allowed the estimation of mprate, prate and ipoint (Fig. 1b) and the second derivative was used to determine the onset of the senescence process (Fig. 1c). The minimum second derivative was calculated for the early and intermediate genotypes but it was not accurately estimated for the late genotypes, therefore it was not considered in this study. However, it can be interpreted as the end of the senescence process and experiments with a longer observation period could make use of this characteristic as well to describe the final stage of senescence.

Repeatability estimation

Once the curve characteristics were calculated, the repeatability of each one was estimated for each year. High values were observed for all of them in the 2 years with values between 0.89 and 0.98. Repeatability was also estimated in each time point and values higher than 0.92 were obtained during the observation period in each year. The high repeatability values observed in the curve characteristics and in the individual time points are due to the experimental design.

G×E interactions

A set of 186 genotypes evaluated in 2004 and 2005 was used to estimate the year effect, the genotype main effect and G×E interaction for prate, mprate, ipoint and onset according to model (3). There was a clear year effect for all four traits (P < 0.001). We observed in 2005 a senescence period lasting longer than in 2004 with a slow progression rate. In 2004, the average mprate and prate were higher (0.36 and 0.86) than in 2005 (0.24 and 0.55), whereas the average onset and ipoint were higher in 2005 (45.97 and 49.90) than in 2004 (37.28 and 39.65). The observation period in 2005 was not long enough for late genotypes to complete the senescence process. Figure 2 shows the senescence fitted curves of three genotypes with different maturity type (CE78: early, CE47: intermediate and CE697: late) in 2004 and 2005 to exemplify the effect of genotype and year in the C×E population. Significant interactions between genotypes and year (P-value < 0.001) were also observed for each curve characteristic. Larger differences in mprate, ipoint and onset were observed between early and late genotypes in 2005, whereas the differences between early and late genotypes were larger for prate in 2004. Early genotypes had a fast senescence development and late genotypes became very late in 2005.

Fig. 2
figure 2

Comparison of senescence development in Finland during 2004 (solid line) and 2005 (dashed line) represented by three genotypes with different maturity type: CE78 (early: red-top line), CE47 (intermediate: blue-middle line) and CE697 (late: yellow-lowest line). The x-axis corresponds to BTT units adjusted for comparison of the 2 years

QTL analysis: characteristics of the curves and single time points

QTLs were detected in each parental map using the individual time points and the characteristics estimated from the curves in each year, considering as a significance threshold a P-value from the Kruskal–Wallis test lower than 0.005. Only the LGs in which QTLs were identified are shown in Fig. 3 and two levels of significance are highlighted (0.001 < P ≤ 0.005, P ≤ 0.001). In the analysis per time point, 10 and 9 time points were used in the QTL analysis in 2004 and 2005, respectively. In the first year, the first time point was excluded from the analysis because all the genotypes were still green. In the analysis using the curve characteristics, although four of them were used for the QTL mapping using KW and IM, results are presented only for onset, prate and ipoint with Kruskal–Wallis. It is due to the fact that high correlations were observed between mprate and prate in the two consecutive years (0.85 and 0.80, respectively) and the same QTLs were detected for both curve characteristics using KW and IM.

Fig. 3
figure 3figure 3

Graphical representation of the QTLs detected in the C (a) and E (b) parental maps. Results from the QTL analysis using data from the field experiments in 2004 and 2005 are presented in the left and right panel respectively. Only the linkage groups in which significant QTLs were detected using individual time points (in 2004: 81,86,91,95,100,105,109,116,123 and 129DAP; in 2005:84,88,93,98,101,106,109,112 and 115 DAP) and curve characteristics (onset, ipoint, prate) are presented. The colors in the linkage groups indicate the significance levels (white: P > 0.005, yellow: 0.001 < P ≤ 0.005, brown: P ≤ 0.001)

Comparing the results using time points and curve characteristics, we observed that QTLs detected in early time points, coincided with QTLs on LG C4, C5A and E3 associated with onset of senescence. QTLs occurring in more advanced stages of the senescence development coincided with QTLs for prate and ipoint as shown in C5 and E5 (they correspond to the same LG in both parents and hence are termed Ch5).

Environment (year)-specific QTLs were found on C1 and C5A associated with prate in 2005 and onset in 2004, respectively. QTLs associated with ipoint were also related with onset in 2005 as it is shown in E7 and Ch5. The pleiotropic QTLs detected for ipoint and onset on both LG explained the high correlation (0.8) observed between the two traits in that year. On the other hand, QTLs consistently found in C4 in the two experimental years were associated with onset and the more significant P-value was observed in 2005 (P < 0.0001).

QTL×QTL interactions

The most significant QTLs for prate, ipoint and onset in each linkage group and all the pairwise interactions between them were used in a linear model to investigate epistatic interactions, within a year. The full model was based on the joined set of main effect QTLs as detected for any of the three curve characteristics. Therefore, the full model was the same for prate, ipoint and onset and it was set according to model (4) including only two-way interactions between QTLs:

$$ \begin{gathered} y = E{\it{3}} + Ch {\it{5}} + E{\it{7}} + C{\it{1}} + C{\it{4}} + C{\it{5}}A + E{\it{3}} \times Ch {\it{5}} + E{\it{3}} \times E{\it{7}} + Ch{\it{5}} \times E{\it{7}} \hfill \\ + C{\it{1}} \times C{\it{4}} + C{\it{1}} \times C{\it{5}}A + C{\it{4}} \times C{\it{5}}A + E{\it{3}} \times C{\it{1}} + E{\it{3}} \times C{\it{4}} + E{\it{3}} \times C{\it{5}}A + Ch{\it{5}} \times C{\it{1}} \hfill \\ + Ch{\it{5}} \times C{\it{4}} + Ch{\it{5}} \times C{\it{5}}A + E{\it{7}} \times C{\it{1}} + E{\it{7}} \times C{\it{4}} + E{\it{7}} \times C{\it{5}}A + {error} \hfill \\ \end{gathered} $$

QTL main effects, QTL×QTL interactions and their corresponding genotypic classes for prate, ipoint and onset are presented in Table 2. To go from the full models to the final models, only significant terms in the full models were retained (P < 0.05). The P-values for the retained terms and adjusted R 2 in the final models for each curve characteristic are presented in Table 3. For instance, after fitting the full model for prate in 2004, 6 terms were retained in the final model (Table 2): two epistatic interactions (E3 × C5a and Ch5 × C4), one main effect QTL (Ch5) and three non-significant main effect QTLs involved in the epistatic interactions (E3, C4 and C5a). In the final model, only the interaction Ch5 × C4 showed to be significant and the main effect of Ch5 and C5a (Table 3).

Table 2 Main effect QTLs and QTL × QTL interactions included in the full models for prate, ipoint and onset of senescence in 2004 and 2005
Table 3 Main effect QTLs and QTL × QTL interactions retained in the final models for prate, ipoint and onset of senescence in 2004 and 2005

In 2004 and 2005, a main effect of Ch5 on prate, ipoint and onset was observed. It is related to a QTL in the same region previously reported as associated with plant maturity type (Celis-Gamboa 2002), which is a trait reflecting the general development of the plant.

In 2005 no epistatic interactions were found between the QTLs but pleiotropic effects of E3 on ipoint and onset, and C4 on prate and onset (Table 3).

In 2004, C4 associated with onset of senescence interacted epistatically with Ch5 associated with prate and ipoint (Table 3). In the interaction C4 × Ch5, eight genotypic classes are present (‘aa ee’, ‘aa ef’, ‘aa eg’, ‘aa fg’, ‘ab ee’, ‘ab ef’, ‘ab eg’, ‘ab fg’) as shown in Table 2 and the average senescence curves per class are presented in Fig. 4. The genotypic classes including ‘ab’ (fitted curves with dashed line) had a faster onset and senescence development than the classes in which ‘aa’ (fitted curves with solid line) was present. The curve ‘ab ee’ showed a higher prate than ‘aa ee’ and the genotypes in this class had in general a faster onset of senescence and died earlier. The classes ‘aa fg’ and ‘ab fg’ were associated with late senescence development and the prate and ipoint of these two classes were the lowest. Genotypes in these two classes were at almost half of the maximum senescence when the observation period ended. Genotypes in the class ‘aa fg’ showed the latest senescence development among all the genotypic classes. Interestingly the class ‘ab fg’ showed the earliest onset but this was followed by a delayed senescence and the genotypes in this group had the latest senescence.

Fig. 4
figure 4

Average senescence curves for each genotypic class in the epistatic interaction between QTLs C4 and Ch5 observed in 2004

Another epistatic interaction was observed in 2004 between C4 associated with onset and E7 related with ipoint (Table 3). Figure 5 shows the average senescence curves observed for the interaction C4 × E7 with eight genotypic classes (‘aa ii’, ‘aa ij’, ‘aa ik’, ‘aa jk’, ‘ab ii’, ‘ab ij’, ‘ab ik’, ‘ab jk’) as shown in Table 2. The classes including ‘ab’ (fitted curves with dotted line) had in general, a faster senescence development than classes in which ‘aa’ (solid line) was present except for the class ‘aa ik’. In this class earlier onset of the senescence process was observed with a faster prate and ipoint while in the group ‘ab ik’ a delayed onset with a slower prate is present. Interestingly, big differences were also observed between the classes ‘aa jk’ and ‘ab jk’. The class ‘aa jk’ showed the same trajectory than the class ‘aa ii’ and both of them have late onset and ipoint making the senescence process slower than in the other classes. On the other hand, the class ‘ab jk’ showed the earliest onset with a delayed senescence that ended up simultaneously with the classes exhibiting the slowest prate (‘aa ii’ and ‘aa jk’).

Fig. 5
figure 5

Average senescence curves for each genotypic class in the epistatic interaction between QTLs C4 and E7 observed in 2004

Discussion

The smoothed GLM procedure in our study allows flexible modeling of different senescence curves in the mapping population. Genotypes with early development showed an s-shape senescence curve, whereas genotypes with slow development exhibited an exponential curve during the observation period. In both cases the curve trajectories would eventually reach the maximum of 7, which is the highest senescence score. In the case of late genotypes, the observation period was not long enough to complete the senescence process. Four characteristics estimated from the curves (onset, prate, mprate and ipoint) allowed us to explain senescence in terms of development. Whereas individual time points only gave us an impression of the trait at selected different moments during the growing season. In the first case sensible biological interpretations could be made and the results were compared with previous studies describing the senescence process in the same population evaluated in the Netherlands (Celis-Gamboa 2002).

The weather and the general conditions of the experiment had a clear effect on the performance of the population in the 2 years of field experimentation in Finland. Clearly the senescence period lasted longer in 2005 and showed a slower progression rate. This was probably due to the higher daily air temperature observed that year. Marinus and Boadlaender (1975) reported that high temperatures accelerate the reproductive and vegetative development of potato plants and delay the senescence process. Our results are therefore in line with their findings.

Comparing the QTL analysis using the individual time points and the curve characteristics, QTLs detected at early time points were associated with onset of senescence. QTLs occurring in intermediate to late time points coincided with QTLs associated with prate and ipoint. Although the results are comparable with both phenotypic data, the analysis using the characteristics of fitted curves offers a biological interpretation of the results in terms of senescence development. In our study, four curve characteristics described the whole process and allowed us to study the progression of the trait. When the C × E population was evaluated in the Netherlands, QTLs on chromosome 5 were identified for different traits related with senescence and life cycle. Malosetti et al. (2006) reported two QTLs on chromosome 5 associated with mid-senescence and rate of senescence. Celis-Gamboa (2002) identified three QTLs on the same chromosome associated with the duration of the plant cycle and one of them was found associated to plant maturity type explaining 40% of the phenotypic variance. In this study, we detected QTLs associated with senescence development on chromosome 5. These are in the same region where the latter QTL was reported by Celis-Gamboa (2002), suggesting a pleiotropic effect of this region for prate, ipoint, onset and plant maturity.

An interesting QTL associated with onset was found on C4, whereas no QTLs have been previously reported for developmental traits on this chromosome. A wide range of functional genes involved in resistance against various pathogens have been reported on chromosome 4, among others a major late blight resistance gene cluster (Celebi-Toprak et al. 2002; Li et al. 1998; Park et al. 2005; Zimnoch-Guzowska et al. 2000).

Epistatic interactions were observed in 2004 between C4 and Ch5 (Fig. 4) and between C4 and E7 (Fig. 5). Interestingly, in both epistatic interactions, the genotypic class ‘ab’ in C4 was associated with early development of the senescence process. In the interaction C4 × Ch5, the genotypic classes ‘ab ee’ and ‘ab ef’ associated with early senescence showed delayed onset but fast prate and ipoint whereas the class ‘ab fg” showed early onset but a senescence period lasting longer. In 2005 no epistatic interactions were observed but main effects of C1 and pleiotropic effects of E3 and C4 on ipoint/onset and prate/onset, respectively.

In summary, by modeling haulm senescence in the C × E potato population using a smoothed GLM procedure, we were able to characterize the curve in terms of developmental traits and to identify QTLs associated with the senescence process. Pleiotropic effects and epistatic interactions between QTLs were detected when two-way interactions were studied. Delayed senescence was associated with particular genotypic classes (C4 × Ch5: ‘ab fg’ and C4 × E7: ‘aa ii’ ‘aa jk’) and some classes showing early onset turned to have delayed senescence development (C4 × Ch5: ‘ab fg’ and C4 × E7: ‘ab jk’).

Different regulatory genes have been reported for onset and progression rate of senescence in Arabidopsis thaliana (Gepstein et al. 2003) and it will be interesting to compare the genes involved in the same process in different crops. The complexity of QTL × E and QTL × QTL × E as well as the performance of genotypes under different environments will be considered in a further study providing insights for the better understanding of adaptation and developmental processes in potato.