A Stochastic Flowering Model Describing an Asynchronically ...

Annals of Botany 90: 405±415, 2002 doi:10.1093/aob/mcf204, available online at www.aob.oupjournals.org

A Stochastic Flowering Model Describing an Asynchronically Flowering Set of Trees F . N O R M A N D 1 , * , R . H A B I B 2 and J . C H A D ê U F 3 BP 180, 97455 Saint Pierre cedex, ReÂunion Island, France, 2INRA, Plante et SysteÁmes de Culture Horticoles, BaÃt. A, Domaine St Paul, Site Agroparc, F-84914 Avignon cedex 9, France and 3INRA, BiomeÂtrie, Domaine St Paul, Site Agroparc, F-84914 Avignon cedex 9, France

1CIRAD-FLHOR,

Received: 11 October 2001 Returned for revision: 18 March 2002 Accepted: 14 June 2002

A general stochastic model is presented that simulates the time course of ¯owering of individual trees and populations, integrating the synchronization of ¯owering both between and within trees. Making some hypotheses, a simpli®ed expression of the model, called the `shoot' model, is proposed, in which the synchronization of ¯owering both between and within trees is characterized by speci®c parameters. Two derived models, the `tree' model and the `population' model, are presented. They neglect the asynchrony of ¯owering, respectively, within trees, and between and within trees. Models were ®tted and tested using data on ¯owering of Psidium cattleianum observed at study sites at elevations of 200, 520 and 890 m in ReÂunion Island. The `shoot' model ®tted the data best and reproduced the strong irregularities in ¯owering shown by empirical data. The asynchrony of ¯owering in P. cattleianum was more pronounced within than between trees. Simulations showed that various ¯owering patterns can be reproduced by the `shoot' model. The use of different levels of organization of the general model is discussed. ã 2002 Annals of Botany Company Key words: Phenology, ¯owering asynchrony, ¯owering model, stochastic model, degree days, Psidium cattleianum.

INTRODUCTION The time course of ¯owering of individual trees and of their population, i.e. the number of ¯owers open daily for the duration of ¯owering, has important consequences for the reproductive success and the genetic structure of a population (Bawa, 1983; Ims, 1990; Murawski and Hamrick, 1991; Hof et al., 1999). In particular, seed and fruit production of self-incompatible species rely on ¯owering times of compatible genotypes overlapping, with important economic consequences for agriculture. The time course of ¯owering of a population can also be considered in relation to the sensitivity of ¯owers to pathogens in order to aid the development of integrated pest and disease management (Dodd et al., 1992), or for human health when the prediction of airborne pollen is used to forecast allergic risks (Frenguelli et al., 1989; Ickovic et al., 1989; Andersen, 1991; Belmonte and Roure, 1991). Plants, in particular tropical trees and shrubs, display a large variety of ¯owering patterns (Gentry, 1974; Bawa, 1983). Asynchronous ¯owering is widespread (Primack, 1980; Bawa, 1983; Ims, 1990; Carthew, 1993), even in mass-¯owering species (Augspurger, 1983). The time course of the population ¯owering, as the superposition of the time course of individual ¯owerings, is affected by the time and duration of individual ¯owering, as well as by the individual number of ¯owers. Two main components can be distinguished in the synchronization of ¯owering (Hof et al., * For correspondence. Fax: +262 2 62 50 58 44, e-mail normand@ cirad.fr

1999): the synchronization between plants and the synchronization within plants. Flowering phenology is mainly under genetic control (Nienstaedt, 1974; Primack, 1980; Mosseler and Papadopol, 1989; Pors and Werner, 1989; Carthew, 1993; O'Brien and Calder, 1993; Boudry et al., 1994; Mitchell-Olds, 1996; Hof et al., 1999), implying that a detailed study of ¯owering phenology must consider plants at the individual level (e.g. Augspurger, 1983; Fripp et al., 1987; Carthew, 1993). Improving synchronization of ¯owering and the ¯owering period have been important objectives for breeding and selection programmes (Janick and Moore, 1975; Hof et al., 1999; Citadin et al., 2001). Flowering phenology is also affected by the environment (Murfet, 1977; Primack, 1980). Temperature is recognized as being the main variable driving the timing of budburst or ¯owering in woody plants, and different models have been proposed to predict the date of ¯owering of a population. They all assume that budburst, or ¯owering, occurs when a critical development threshold a is reached, the stage of development being a sum of daily rates of development. The models differ in the expression of the rate of development as a function of temperature (e.g. HaÈnninen, 1987; Chuine et al., 1998, 1999), and in the way they integrate chilling requirements for temperate trees (Cannell and Smith, 1983; HaÈnninen, 1987; Murray et al., 1989; Kramer, 1994; Chuine et al., 1998, 1999). The predicted date of ¯owering is the day on which the critical development threshold is reached. It generally corresponds to the onset of bloom, to the mid-bloom date, or to the date of maximum concentration of airborne pollen (Boyer, 1973; Chuine et al., 1998, 1999). ã 2002 Annals of Botany Company

406

Normand et al. Ð Modelling Asynchronous Flowering

However, very little attention has been paid to modelling of the time course of ¯owering of a tree or a population. Agostini et al. (1999) proposed a stochastic ¯owering model for an orchard of kiwifruit [Actinidia deliciosia (A. Chev.) C.F. Liang & A.R. Ferguson] female vines (see Materials and Methods). In such a clonal population, the time course of individual ¯owering is expected to be reasonably similar, with a strong synchronization of ¯owering between plants. If ¯owering is a regular process on individuals and the environment is homogeneous, then the time course of ¯owering for the population and for individuals should not differ, except for the number of ¯owers open per day. More generally, in a non-clonal population, such a model would be inappropriate due to asynchronous ¯owering and differences in the time course of individual ¯owerings (e.g. Primack, 1980; Augspurger, 1983; Fripp et al., 1987; Carthew, 1993). Simple empirical functions have been proposed to describe the time course of ¯owering at the individual level using plant-speci®c parameters (Fripp et al., 1987; Medan and Bartoloni, 1998; Hof et al., 1999), but without integration at the population level. With regard to classical phenological models, it appears necessary to deal with ¯owering variability both within and between trees to model the time course of ¯owering at the tree and at the population level. The within-tree ¯owering variability, i.e. the source of the time course of ¯owering at the tree level, may be related to genetic factors, to physiological factors, to the time-lag in ¯ower induction or in bud break, or to the buds' effective temperature. The between-trees variability is related to genetic and local environmental factors. The study and modelling of ¯owering phenological processes is a way to build an explanatory model to simulate the time course of ¯owering. However, results may be speci®c to the species studied and such models will probably need inputs that are dif®cult to obtain and will therefore be of limited use, in particular for prediction. Moreover, data on these processes are lacking. Another approach, that we chose to adopt, is to model the components of ¯owering variability in order to take them into account explicitly in the ¯owering model. The objective was not to explain the processes underlying the time course of ¯owering, but to simulate their effects with a model of general use. This paper proposes a theoretical model to simulate the time course of ¯owering of individual trees and their population which integrates the ¯owering variability between and within trees. A general model is presented and hypotheses are proposed to simplify it and reduce the number of parameters. Derived models neglecting the within-tree or the within- and between-trees variability are also presented. A procedure to estimate the model parameters is proposed. The models are then ®tted and tested using experimental data obtained in the subtropical ReÂunion Island on strawberry guava (Psidium cattleianum Sabine) populations growing in natural conditions at three elevations. The simplifying assumptions used for the model construction are tested on P. cattleianum. Model structure and uses are then discussed.

MATERIALS AND METHODS Model description

The heart of the model was developed by Osawa et al. (1983) to describe bud phenology of balsam ®r [Abies balsamea (L.) Mill.]. Dennis et al. (1986) adapted it to insect development, and Agostini et al. (1999) to kiwifruit ¯owering. We recall here its assumptions and formulation. It assumes that the development of a given bud, here a ¯ower bud, is a stochastic process consisting of the accumulation of small increments of development ud. This process is supposed to begin with the ¯ower bud at a given stage at time t = 0. The process S(t) is de®ned as the amount of development time a ¯ower bud has accumulated by actual time t: X ud : St dt

The amount of development time, S(t), and t are expressed in degree days. If the increment of time Dd is small and the increments of development ud are independent and identically distributed with an expectation E(ud) = Dd, then S(t) is normally distributed with mean t and variance s2t (Osawa et al., 1983). Flowering occurs when the amount of development S(t) passes through a threshold a. So the probability that a given ¯ower bud is ¯owering at time t is: aÿt 1 pt ProbSt > a Y p t where 1

Yu u

1 x2 p eÿ 2 dx 2

is the cumulated probability between u and +` of the standard normal distribution. The probability that a given ¯ower bud is ¯owering between ti ± 1 and ti is: a ÿ ti a ÿ tiÿ1 2 Pti Y p ÿ Y p ti tiÿ1 The implicit assumption of this model is that all the ¯ower buds follow the same stochastic process, i.e. they accumulate increments of development in the same way and have the same ¯owering threshold. Let us now consider a population of K trees, and a development time scale calculated at the population level. The ¯owering variability between trees is related to different individual ¯owering patterns linked to the genetic variability and also, with reference to a common development time scale, to the differences of effective temperature at the level of each tree. We assume that the ¯owering variability within trees is related to different ¯owering patterns of the tree ¯owering units. We de®ne the ¯owering unit as a level of organization whose ¯ower buds are in a

Normand et al. Ð Modelling Asynchronous Flowering similar stage of development and experience a similar environment. Consequently, they have a similar pattern of development. For convenience, we will later call the ¯owering units ¯owering shoots. However, the architectural level of a ¯owering unit needs to be determined for each species (e.g. in¯orescence for some species, etc.). The within-tree ¯owering variability is de®ned by two components: (1) the ¯owering variability of the ¯ower buds on a ¯owering unit, described by the stochastic process [eqn (2)]; and (2) the variability in ¯owering time between ¯owering units. The ¯owering variability can be formulated as a variable ¯owering threshold a and a variable parameter of variance relative to ¯owering s2 at the shoot and at the tree levels. Increasing values of the parameter of variance relative to ¯owering are expected when the level of organization increases (shoot < tree < population). In the population of K trees, each tree, p, has FSp ¯owering shoots. A ¯owering shoot, denoted fsj,p j 2 1; FSp , bears Nj,p ¯ower buds that are supposed to follow the same stochastic process of development. Equation (2) can then be applied at the ¯owering shoot level with the speci®c parameters aj,p and sj,p2 common to all its ¯ower buds. Let us consider the expectation of parameters aj,p, ap = E(aj,p) j 2 1; FSp , as the tree ¯owering threshold for tree p. The difference between the ¯owering threshold of the ¯owering shoot and the tree ¯owering threshold, dj,p = aj,p ± ap, expressed in thermal time units, represents the advance (dj,p < 0) or the delay (dj,p > 0) of the mid-bloom on the shoot fsj,p compared with the mid-bloom on all the shoots of the tree. By construction, the expectation of dj,p is zero. Therefore, the probability that a ¯ower bud on shoot fsj,p, j 2 1; FSp ; p 2 1; K, will ¯ower between ti ± 1 and ti is:

and extend its use, in particular for prediction, we propose to model the variability of these parameters using statistical distributions whose parameters are estimated from the data. The components of the ¯owering variability are thus quanti®ed. We propose three hypotheses for the distributions of ap, dj,p and sj,p: (1) the differences dj,p follow the same normal distribution N(0,d2) for all trees, i.e. the tree has no in¯uence on the distribution of shoot ¯owering thresholds around the tree ¯owering threshold; (2) the parameter of variance relative to ¯owering at the shoot level, sj,p2, is independent of the shoot and the tree, and is called ss2; and (3) the tree ¯owering thresholds, ap, follow a normal distribution N(as,bs2). Given the ®rst and second hypotheses, we can calculate the probability of ¯owering between ti ± 1 and ti for a ¯ower bud on tree p, unconditional on the shoot, using eqn (3): ap u ÿ ti p ÿ Pp ti Y s ti ÿ1 ap u ÿ tiÿ1 p gudu hap ; ti Y s tiÿ1 1

K 2

pX K

FSp

p1

depends on the number of shoots and trees. This general formulation is not easy to use for prediction, and parameter values are speci®c to the shoots and trees upon which they are estimated. These parameters are the values of the tree ¯owering threshold ap for each tree and of the two parameters of the stochastic ¯owering process (aj,p = ap + dj,p and sj,p) for each ¯owering shoot. To simplify the model

4

where u2 1 gu p eÿ 22 2

is the density function of the normal distribution N(0,d2). Given the third hypothesis, we can calculate the probability of a ¯ower bud ¯owering between ti ± 1 and ti, unconditional on the shoot and the tree, using eqn (4): 1

ap dj; p ÿ ti ap dj; p ÿ tiÿ1 p p ÿY 3 Pj; p ti Y j; p ti j; p tiÿ1 This is the general form of the ¯owering model, where ap, dj,p and sj,p2 are parameters. All ¯ower buds on shoot fsj,p have the same probability of ¯owering between ti ± 1 and ti. This probability depends only on time and on the model parameters. Therefore, the model assumes that ¯owering of a given ¯ower bud on shoot fsj,p is not affected by the other ¯ower buds on the shoot, i.e. the ¯ower buds are independent with regard to their probability of ¯owering. The model parameters are numerous and their number

407

Pti

1

Pp ti f ada ÿ1

ha; ti f ada

5

ÿ1

where f a

aÿs 2 1 ÿ p e 2 s 2 s 2

is the density function of the normal distribution N(as,bs2). Or: a u ÿ ti p ÿ Y Pti s ti ÿ1 ÿ1 a u ÿ tiÿ1 p gudu f ada Y s tiÿ1 1 1

6

where g(u) and f(a) are as in eqns (4) and (5), respectively. This formulation of the model is of general use and has only four parameters related to the components of the ¯owering variability: ss2 is the parameter of variance relative to ¯owering at the shoot level and expresses the time course of the ¯owering process at the shoot level; d2 is

408


the variance of the central normal distribution of the differences dj,p on the tree and expresses the synchronization of shoots ¯owering within a tree (the higher the value of d2, the more asynchronous is ¯owering among shoots on a tree); as is the mean of the normal distribution of the tree ¯owering thresholds, i.e. the mean ¯owering threshold of the ¯owering shoot population taking into account the ¯owering of individual shoots; and bs2 is the variance of the normal distribution of the tree ¯owering thresholds and expresses the ¯owering synchronization between trees (the lower the value of bs2, the more synchronous is ¯owering between trees). To evaluate the relevance of this model, called hereafter the `shoot' model, two simpler models have been derived from it for comparison. The ®rst, called the `tree' model, considers that the within-tree variability can be neglected (i.e. the ¯owering units of a tree are quasi-synchronous): dj,p = 0. The effect of trees on ¯owering is then characterized by their respective tree ¯owering threshold ap¢, distributed according to a normal distribution N(at, bt2). We hypothesize that the parameter of variance relative to ¯owering, st2, is independent of trees. This model has three parameters: at, bt2 and st2, and the probability of a ¯ower bud ¯owering between ti ± 1 and ti, unconditional on the tree, is: 1

Pti ÿ1

a ÿ ti a ÿ tiÿ1 Y p ÿ Y p f ada t ti t tiÿ1

7

where f a

aÿt 2 1 ÿ p e 2 t 2 t 2

is the density function of the normal distribution N(at, bt2). The second simpler model, termed the `population' model, considers that there is no effect of the trees on ¯owering, i.e. that their respective ¯owering is quasisynchronous at the population level. Trees have the same ¯owering threshold a, and there is a parameter of variance relative to the ¯owering of the whole population: s2. This model, similar to that of Agostini et al. (1999), has two parameters, a and s2, and is described by eqn (2). Model simulation

The number of ¯ower buds on each ¯owering shoot must be known for simulations using the `shoot' model, so that their respective weight in the total ¯owering can be taken into account. A development time scale is also required. Simulations are realized as follows. First, a ¯owering threshold, ap, is randomly sampled for each tree in the normal distribution N(as,bs2). Then, a difference, dj,p, is randomly sampled for each ¯owering shoot of each tree in the normal distribution N(0,d2). Given ap, the ¯owering threshold aj,p of each ¯owering shoot is calculated: aj,p = ap + dj,p. Given aj,p, the ¯owering probability distribution (FPD) of the Nj,p ¯ower buds of the shoot is determined

using eqn (3), using sj,p or a common ss according to the model used. A ¯owering date is randomly sampled within FPD for each of the Nj,p ¯ower buds of the shoot. Simulated ¯owering dates are then aggregated at the shoot, tree or population level. To illustrate the effect of between- and within-tree synchronization on the ¯owering pattern, ¯owering of a ®ve tree population has been simulated with four distinct cases of synchronization: A, asynchrony strong between trees and weak within trees (bs = 250 °Cd, d = 2 °Cd); B, strong asynchrony between and within trees (bs = 150 °Cd, d = 150 °Cd); C, asynchrony weak between trees and medium within trees (bs = 20 °Cd, d = 100 °Cd); D, weak asynchrony between and within trees (bs = 20 °Cd, d = 2 °Cd). The other parameters were set at as = 1000 °Cd and ss = 0´7 °Cd. The development time scale ran from 500 to 1500 °Cd, with a 10 °Cd daily increment. Each tree had the same weight: 100 ¯owering shoots, each with 15 ¯ower buds, i.e. 1500 ¯ower buds per tree. Data

Flowering of Psidium cattleianum (strawberry guava tree) was recorded from November 1998 to January 1999 at three sites at elevations of 200, 520 and 890 m on the windward east coast of the subtropical ReÂunion Island, Indian Ocean (21°06¢S, 55°32¢E). The sites at elevations of 200 and 890 m were fallow lands invaded by feral strawberry guava trees. The trees at 520 m were from an experimental 5-year-old seedling orchard which had not received any cultural care for the previous year and were thus considered feral. All individuals were genetically different on all the sites. On each site, ten trees of similar size (1´5±2 m tall) were sampled. Flower buds are borne on new shoots emerging from the terminal branches at the end of the cool, dry season. Generally, one to three shoots emerge quasisimultaneously from the same terminal branch (Normand and Habib, 2001). These shoots and their ¯ower buds are therefore the same age and experience the same microenvironment. Consequently, we considered all the ¯owering shoots emerging from a terminal branch as the ¯owering unit for the purposes of the model, termed hereafter the `¯owering shoot'. Terminal branches were sampled randomly in the tree canopy before bud burst: 20 branches were sampled on eight trees and 100 on two trees at the site at 520 m and 20 per tree at the sites at 200 and 890 m. On each tree on which 100 shoots were observed, 20 shoots were randomly sampled among the 100 shoots, and parameters were estimated using these 20 sampled shoots to avoid a higher weight of these trees in the likelihood estimation due to more ¯owers being observed. On the other hand, the complete 100-shoot data sets were used to estimate the distribution of parameters dj,p and sj,p on a large sample to test the simplifying hypotheses made to derive the `shoot' model. Each ¯ower stays open for 1 d: the ¯ower bud bursts in the morning (0700±0900), and anthers and petals dry and fall in the late afternoon. At each site, open ¯owers were counted daily on the shoots that had emerged from the terminal branches sampled. The data set obtained from the

Normand et al. Ð Modelling Asynchronous Flowering site at 520 m was used for parameter estimation and model comparison, and the data sets from the sites at 200 and 890 m were used for model validation. Mean daily air temperatures, (Tmax + Tmin)/2, were estimated at the sites at 520 and 890 m using the daily thermic gradient between temperatures recorded at 40 and 1025 m, as described by Normand and Habib (2001). Mean daily air temperatures were determined directly at the site at 200 m. The development time scale used was a thermal time scale (Cannell and Smith, 1983). Data collected in orchards over several years and at several elevations showed that the delay between the triggering of a ¯owering shoot upon fertilization and mid-bloom was 914 °Cd with a base temperature, Tb, of 8´1 °C. Although ¯owering data for this study were collected on natural, non-triggered ¯owers, heat units were calculated using 8´1 °C as the base temperature, and the starting date, t0, for the sum of heat units was determined at each site so that mid-bloom occurred on the day when calculated degree days were the closest to 914 °Cd. Empirical and simulated time courses of ¯owering were expressed as ¯owering frequency distribution (FFD). The ¯owering frequency on day ti was: Fti

nti N

where n(ti) was the number of ¯owers open on day ti, and N was the total number of ¯owers observed. Statistical analysis

For each model, parameters were estimated by the maximum likelihood method (DacuÈnha-Castelle and Du¯o, 1982) on data giving the time course of the number of open ¯owers on a sample of ¯owering shoots observed on a sample of M trees. Each tree, p, has FSOp ¯owering shoots, denoted fsoj,p, each bearing NOj,p ¯ower buds j 2 1; FSOp ; p 2 1; M. Observations were made from t1 to tD and covered the whole ¯owering period of the M trees. At time ti i 2 1; D, nj,p(ti), ¯ower buds burst into bloom on shoot fsoj,p. Equation (3) assumes that ¯ower buds are independent with regard to their probability of ¯owering, so that the probability of the observed time course of ¯owering from t1 to tD on shoot fsoj,p follows a multinomial distribution and is D NOj; p ! iY Pj; p ti nj; p ti iQ D nj; p ti ! i 1 i1

The multinomial coef®cient depends only on the sample size and does not in¯uence the parameter estimation and the maximum likelihood tests. It is then removed from the following calculations. Furthermore, we hypothesize that ¯owering shoots on a tree are independent with regard to their probability of ¯owering. Then, the likelihood of the sample for the shoot model is given by:

L

pY K p1

409

1 j FSO 1

Y

ÿ1

j1

a u ÿ tiÿ1 p Y s tiÿ1

p

ÿ1

a u ÿ ti p ÿ Y s t i

iY D i1

nj; p ti

gudu f ada

8

where 1

Yx x

and

y2 1 p eÿ 2 dy; 2

f a

u2 1 gu p eÿ 22 2

a ÿ s 2 1 ÿ p e 2 s 2 s 2

For the `tree' and the `population' models, the likelihood of a sample is a simpler expression. The likelihood estimates of the parameters ss2, d2, as and 2 bs are the values which maximize L, or minimize ±log(L). Calculations were performed with a minimization program for non-linear functions (Splus 2000 statistical package) using a general quasi-Newton optimizer (Mathsoft, 1999) on the M = 10 trees (2370 ¯ower buds) observed at the site at 520 m. Computing limitations did not permit the minimization of ±log(L) for the `shoot' model including two nested integrals [eqn (8)]. We estimated the parameters (as, bs, d, ss) by selecting the value that minimized ±log(L) in a range of values centred on a ®rst estimation derived from the data (as, d and ss were the means of these parameters estimated at the tree level, and bs was the standard deviation of the tree ¯owering thresholds). Although not optimal, this method leads to asymptotically unbiased parameters whose variance tends to zero when the number of ¯owers tends to in®nity. The maximum likelihood test (Lehman, 1983) was used to compare the three models of increasing complexity (`population', `tree' and `shoot'), and to test the pertinence of some of the hypotheses made in model construction. The test evaluates the increase in likelihood between two models compared with the difference in their numbers of parameters. For a pair comparison of n models, twice the difference of log-likelihood, calculated with the estimated parameters of a data set, was compared with a c2 value for P = 5/n % and a number of degrees of freedom equal to the difference between their respective parameters numbers. Distributions were compared using a bilateral Kolmogorov± Smirnov's test (Sprent, 1989).

RESULTS Model simulations

Simulations performed with the `shoot' model in Fig. 1 show four levels of between- and within-tree synchronization. The main ¯owering patterns were well repre-

410


F I G . 1. Patterns of ¯owering, simulated using the `shoot' model, of a set of ®ve trees each bearing 100 shoots with 15 ¯ower buds per shoot for different levels of ¯owering synchronization between (bs) and within (d) trees: A, bs = 250 °Cd, d = 2 °Cd; B, bs = 150 °Cd, d = 150 °Cd; C, bs = 20 °Cd, d = 100 °Cd; and D, bs = 20 °Cd, d = 2 °Cd. The other parameters are set at as = 1000 °Cd, ss = 0´7 °Cd.

sented by the model (Gentry, 1974; Bawa, 1983), in particular, the case of a strong asynchrony between trees, which is not reproduced by classical ¯owering models (Fig. 1A). Moreover, these patterns may be modi®ed by differing numbers of ¯ower buds on the shoots and trees. Flowering of individual trees was also simulated (data not shown), and can be useful to quantify the degree of overlap in ¯owering among trees (Primack, 1980; Augspurger, 1983). Model ®tting and model comparison

The data set collected at the site at 520 m is summarized in Table 1. Not all the terminal branches studied produced ¯owering shoots. The number of ¯ower buds per ¯owering shoot varied from one to 50. The total number of ¯ower buds observed was 4151, and 2370 ¯ower buds were taken into account for parameter estimation following a 20-shoot sampling among the 100 shoots of two trees (Table 1). At the sites at 200 and 890 m, 695 and 972 ¯ower buds were observed. Differences were related to different ¯ushing and ¯owering intensities of the trees among the sites. Only eight trees ¯owered at the site at 200 m. Parameters of the three models ®tted on the 520 m data set were: `shoot' model as = 927´3 °Cd, bs = 40´6 °Cd, d = 95´4 °Cd, ss = 0´96 °Cd; `tree' model at = 923´7 °Cd, bt = 50´8 °Cd, st = 2´61 °Cd; `population' model a = 914´4 °Cd, s = 3´12 °Cd. As expected, the value of the parameter of variance relative to ¯owering increased with the level of organization of the model (shoot < tree < population). It re¯ected the

TA B L E 1. Numbers of ¯owering shoots and ¯ower buds, and the onset and end of ¯owering of the sample of ten Psidium cattleianum trees (1±10) at the site at 520 m in ReÂunion Island Number of Number of Onset of End of Tree ¯owering shoots ¯ower buds ¯owering (°Cd) ¯owering (°Cd) 1 2 3 4 5 6 7 8 9 10 3s* 6s*

18 18 92 19 19 94 20 17 20 13 20 17

287 185 1112 228 275 1071 314 286 301 92 233 169

621´4 881´3 710´2 821´3 721´0 731´8 710´2 857´6 721´0 832´9 721´0 743´1

977´7 1153´2 1203´2 1140´6 1216´2 1190´3 1039´5 1255´6 1270´3 1255´6 1078´3 1190´3

* 3s and 6s correspond to a 20-shoot random sampling among the 100 tagged shoots of trees 3 and 6, respectively.

dispersion of ¯owering at each level. bt and bs were high, indicating the asynchrony of ¯owering between trees (Table 1). Individual tree ¯owering thresholds, ap¢, estimated by the `tree' model, varied from 837´7 to 998´9 °Cd. The difference represented 14 d. Likewise, d was high, indicating the asynchrony between shoots on each tree. The largest difference between shoot ¯owering thresholds, aj,p, on a tree was 445 °Cd, or 38 d.


411

TA B L E 2. Pair comparisons of the `shoot', `tree' and `population' ¯owering models by tests of maximum likelihood, each at P = 1´7 %, with the ¯owering data of ten Psidium cattleianum trees from a site at 520 m in ReÂunion Island. Model 1

Model 2

Shoot Shoot Tree

Tree Population Population

2DlL

d.f.

c20.017

P

3896´8 4695´0 798´2

1 2 1

5´7 8´1 5´7

0 0 0

2DlL, Twice the difference of log-likelihood of the models; d.f., difference between the number of parameters of the models; c20.017, test's threshold (0´983 quantile of the cumulated c2 distribution with d.f. degrees of freedom); P, 1 ± cumulated probability of 2DlL on a c2 distribution with d.f. degrees of freedom.

The sample log-likelihood was high for each model because of the large sample size. The test of maximum likelihood was highly signi®cant for all pair comparisons of the three models (Table 2). The large gain of likelihood between the `shoot' model and the two other models indicated that the former ®tted the data in a more likely manner than the other models, i.e. taking into account the shoot level brought a large gain of accuracy in describing the data. To study the behaviour of each model at the population level, 1000 simulations were run, using shoot and/or tree ¯owering thresholds randomly sampled in estimated distributions, and using the thermal time scale and the number of ¯ower buds per shoot and tree observed at the site at 520 m. The 0´025 and 0´975 quantiles and the median of the simulated ¯owering frequencies were determined for each date, ti. The empirical data were then graphically compared with the range of 95 % of the simulated values belted by these quantiles (Fig. 2). Empirical FFD showed strong irregularities, particularly during peak ¯owering, indicating that ¯owering was not a smooth process. These irregularities were due to the simultaneous ¯owering of several shoots bearing many ¯ower buds. The three models gave a good general ®t of the ¯owering pattern, in particular the midbloom date and the overall ¯owering duration. The median FFDs simulated by each model were not signi®cantly different (Kolmogorov±Smirnov's test for pair comparison, P > 0´84), indicating that, on average, the three models gave similar FFDs (data not shown). The models differed in the width of the inter-quantiles band, which traduces the capacity of the model to simulate irregularities in the time course of ¯owering as those of the empirical FFD. The band was narrower as the level of organization of the model increased (Fig. 2). Consequently, the empirical FFD was more satisfactorily included in the `shoot' and `tree' model inter-quantile bands than in that of the population model. The `tree' model band-width was as large as that of the `shoot' model during peak ¯owering, but was narrower at the onset and the end of ¯owering. The `tree' model was thus as ef®cient as the `shoot' model at representing the strong irregularities of empirical FFDs during peak ¯owering, but was less ef®cient at the onset and the end of ¯owering. This

F I G . 2. Empirical time course of ¯owering of ten P. cattleianum trees (2370 ¯ower buds) at the site at 520 m elevation in ReÂunion Island (solid line), and 0´025 and 0´975 quantiles of 1000 simulations (dashed lines) run with the `shoot' model (A), the `tree' model (B) and the `population' model (C) ®tted on the empirical data.

was probably the consequence of the random attribution, by the `shoot' model, of a ¯owering threshold to each ¯owering unit whose weight in the tree ¯owering is its number of ¯ower buds. Thus, ¯owering units with numerous ¯ower buds could ¯ower early or late, and increased the simulated ¯owering frequency at the onset or at the end of overall ¯owering.

412

Normand et al. Ð Modelling Asynchronous Flowering TA B L E 3. Kolmogorov±Smirnov's tests for the central normal distribution of the differences d j,p between the shoot ¯owering thresholds and the tree ¯owering threshold of ten Psidium cattleianum trees at the site at 520 m in ReÂunion Island

F I G . 3. Validation of the `shoot' model at 200 m of elevation: empirical ¯owering frequency distribution (solid line) with the 0´025 and 0´975 quantiles of 1000 simulations (dashed lines) run with the `shoot' model ®tted at 520 m. The 200 m thermal time scale, and the empirical numbers of ¯ower buds per shoot were used for the simulations.

F I G . 4. Validation of the `shoot' model at 890 m of elevation: empirical ¯owering frequency distribution (solid line) with the 0´025 and 0´975 quantiles of 1000 simulations (dashed lines) run with the `shoot' model ®tted at 520 m. The 890 m thermal time scale, and the empirical numbers of ¯ower buds per shoot were used for the simulations.

Model validation

When expressed in thermal time, the duration of overall empirical ¯owering at the three sites was similar (592´8 °Cd at 200 m, 648´9 °Cd at 520 m and 472´2 °Cd at 890 m), indicating that ¯owering was mainly driven by temperature, which was our basic hypothesis. However, ¯owering patterns were different, with a straight ¯owering peak during mid-bloom at the site at 200 m (Fig. 3), and just before mid-bloom at the 890 m site (Fig. 4). These peaks were related to heavily ¯owering trees with similar ¯owering thresholds. At the 200 m site, four trees (of eight) bore 71´1 % of the observed ¯ower buds and had very similar ¯owering thresholds (902´7±918´7 °Cd). At the site at 890 m, three trees (of ten) bore 56´8 % of the observed ¯ower buds and also had similar ¯owering thresholds (861´6±933´4 °Cd).

Tree

ks

P

Estimated dp (°Cd)

Calculated dp (°Cd)

1 2 3 4 5 6 7 8 9 10 3s* 6s*

0´1867 0´1300 0´0962 0´1194 0´0991 0´0501 0´1230 0´1385 0´1630 0´1210 0´0891 0´1488

0´50 0´88 0´36 0´92 0´98 0´97 0´89 0´86 0´61 0´98 0´99 0´79

64´0 55´1 67´8 81´0 91´8 86´1 70´3 70´0 110´3 113´8 75´1 95´1

66´0 58´4 68´6 83´8 95´5 87´1 72´4 73´9 114´4 118´6 77´6 98´6

The standard deviation (dp) of the central normal distributions, estimated by the model or calculated from the data, are presented. Number of ¯owering shoots per tree are as in Table 1. ks, Kolmogorov±Smirnov's test threshold; P, test's probability. * 3s and 6s correspond to a 20-shoot random sampling among the 100 tagged shoots of trees 3 and 6, respectively.

The `shoot' model was validated for temperature and site variation at the 200 and 890 m sites using the previous method. Simulations were run using parameter values estimated from the data set collected at the 520 m site, and the empirical thermal time scale and number of ¯ower buds per shoot relating to each site (Figs 3 and 4). Empirical FFDs were satisfactorily included in the range of 95 % of the simulated values, except for the ¯owering peaks observed at these sites. The model gave a good representation of overall ¯owering duration and mid-bloom date at both sites. The `tree' and `population' models were validated in the same way. At each site, results of model behaviour were similar to those of the model comparison at 520 m (data not shown). Hypothesis testing

Three hypotheses were made to simplify the general form of the ¯owering model given by eqn (3) and to derive the shoot model [eqn (6)]. These hypotheses were tested on the 520 m data set to verify their biological relevance for P. cattleianum. Kolmogorov±Smirnov's tests indicated that the differences, dj,p, followed a central normal distribution for all trees, but with different standard deviations, dp (Table 3). The hypothesis of the same normal distribution of the differences dj,p on all trees was tested using a maximum likelihood test comparing a model (ap¢,sp2,d) with the same d for all trees with a model (ap¢,sp2,dp) with dp estimated for each tree. Each tree was characterized by the estimated values of its ¯owering threshold, ap¢, and its parameter of variance relative to the ¯owering of its shoots sp2. The test was not signi®cant at the 5 % level (2DlL = 16´6, d.f. = 9, c2 = 16´9, P = 0´06), indicating that this hypothesis was veri®ed

Normand et al. Ð Modelling Asynchronous Flowering by our sample. However, the test's probability was closed to the signi®cance threshold, suggesting that further work is necessary to verify, and eventually modify, this hypothesis for P. cattleianum. The parameter of variance relative to ¯owering at the shoot level, sj,p2 [eqn (3)], was hypothesized to have a unique value, ss2, independently of the shoot and the tree. We ®rst tested its independence with the shoot using a test of maximum likelihood comparing a model (aj,p, sj,p2) in which each shoot was characterized by the estimated value of its ¯owering threshold aj,p and its parameter of variance relative to ¯owering sj,p2, with a model (aj,p, sp2) in which the shoots of each tree had a common parameter of variance sp2 relative to ¯owering. The test was highly signi®cant (2DlL = 1142, d.f. = 71, c2 = 202´5, P = 0´00). The hypothesis of the same sp2 for all shoots of a tree was therefore rejected, and a fortiori the hypothesis of the same parameter of variance relative to ¯owering ss2 for all shoots of the tree population. The sj,p2 values were highly variable among trees (sj,p coef®cient of variation on trees 3 and 6: 72´2 and 74´2 %, respectively). The sj,p2 traduced, in part, the variability in development time accumulation among the ¯ower buds of a ¯owering shoot. Two sources of variation could be put forward. First, slight signi®cant positive correlations were found between sj,p and the number of ¯owers of the shoot (r = 0´37, n = 89, P < 0´001 for tree 3, and r = 0´52, n = 93, P < 0´001 for tree 6). This suggested a differential distribution of assimilates among ¯ower buds when they are numerous on a shoot, leading to variable rates of development, and a larger sj,p on the shoot. Secondly, variable sj,p2 among shoots of a tree could be related to differences in the effective temperature of their ¯ower buds, linked to the shoot position in the canopy (inside/outside, cardinal orientation). The ¯owering thresholds ap¢ of the trees at the 520 m site estimated by the tree model followed a normal distribution with mean at and variance bt2 (Kolmogorov±Smirnov's test, n = 10, ks = 0´2129, P = 0´68). DIS CUS S ION Model structure

A stochastic model simulating the time course of ¯owering of an asynchronically ¯owering population of trees has been presented. It is based on a synthetic representation of the ¯owering variability. Flowering at a particular level (tree or population) is the aggregation of ¯owering of the lowerlevel components (¯owering unit or tree, respectively). The basic component is the ¯owering unit, whose ¯ower buds are supposed to follow the same stochastic process of development. The ¯owering unit must be determined for each species. The model considers the ¯owering variability both between and within trees. Within-tree variability is related to the time course of ¯owering at the ¯owering unit level and to differences in ¯owering time between ¯owering units. Between-tree variability is related to different tree ¯owering thresholds due to genetic and environmental factors. The inputs are a development time scale and the number of ¯ower buds per ¯owering unit or per tree,

413

obtained by observation or as plant growth model outputs. The aggregative structure of the model allows different levels of detail in the outputs, from the ¯owering unit to the population time course of ¯owering. The general form of the model [eqn (3)] has a large number of parameters speci®c to the shoots and trees upon which they are estimated and is not relevant for use on other trees. Three hypotheses on the distribution of the parameter values were proposed to simplify the model and derive the `shoot' model, which requires only four parameters quantifying the components of the ¯owering variability. These hypotheses have biological meaning and their pertinence should be tested for each species. The in¯uence of environmental and genetic factors on ¯owering phenology tends to be under-estimated by the hypotheses as the individual tree has no in¯uence on the distribution of the differences dj,p and on the parameter of variance relative to ¯owering at the shoot level. Each tree is characterized by its own ¯owering threshold and by the number of ¯ower buds on its shoots. Nevertheless, the simulated time course of ¯owering of n trees is not the superposition at different times of n ¯owering of the same tree. The model randomly samples, in an estimated distribution, a ¯owering threshold for each ¯owering shoot whose weight is its number of ¯ower buds. It then reproduces the differences in ¯owering phenology among trees and among shoots within a tree, in the limits of the variability of the data used to estimate the model parameters. These data must therefore be representative of the existing variability, in particular with respect to the in¯uence of genetic and environmental factors, in order to use the model on other populations. The simulations performed with different levels of ¯owering synchronization between and within trees show that the main ¯owering patterns are well represented by the `shoot' model (Fig. 1). For example, Gentry (1974) found ®ve distinct ¯owering patterns among the American Bignoniaceae, mainly related to different strategies with regard to pollinators (see also Bawa, 1983; Ims, 1990). They correspond to various levels of synchronization within and between trees, and can be well simulated by the model: Fig. 1B corresponds to the `cornucopia' ¯owering pattern and Fig. 1D to the `big band' ¯owering pattern (Gentry, 1974). Two simpler models were derived from the `shoot' model: the `tree' model that neglects within-tree ¯owering asynchrony, and the `population' model that also neglects the between-tree ¯owering asynchrony (e.g. Agostini et al., 1999) and considers that all the ¯ower buds of the population follow the same stochastic process of development (i.e. no genetic or environmental in¯uence). Model adaptability

The model can be adapted to the biological characteristics of a particular species or environment, and to the objectives and needs of a particular study. The hypotheses made to derive the `shoot' model may not be suited for some species. It is then possible to use different distributions for the differences dj,p and the tree ¯owering thresholds in eqns (4) and (5), or to ®x a tree ¯owering sequence in the case of a strong environmental (Primack, 1980) or genetic (Primack,

414


1980; Pors and Werner, 1989; Carthew, 1993; O'Brien and Calder, 1993; Mitchell-Olds, 1996; Hof et al., 1999) in¯uence on ¯owering phenology. Likewise, a particular distribution of the parameter of variance relative to ¯owering at the shoot level can be used. The analytic formulation of the `shoot' model would then be more complex. The model can be used to simulate the time course of ¯owering for a predictive or descriptive purpose, but the accuracy of the simulation is not the same. The most accurate way to reproduce the observed time course of ¯owering is to estimate the parameters at the tree and/or ¯owering unit level and to ®x their value in the simulation process. The variability is then only induced by the stochastic process of ¯ower bud development. If the objective is to predict the time course of ¯owering, then simulations are run with the ¯owering thresholds randomly sampled in distributions whose parameters are also estimated. Simulations are less accurate, but the model is of more general use, provided that the environment and the population whose time course of ¯owering is simulated are included in the domain of validity of the model de®ned by the data upon which parameters were ®tted. Our study on P. cattleianum showed that overall ¯owering duration and mid-bloom date are well estimated, even if parameters are ®tted in another place. Moreover, the range of the most likely simulated data can be estimated from a large number of simulations (Figs 2±4). A limitation of the predictive ability of the model may arise from a strong lack of balance in the number of ¯ower buds of the components of a ¯owering level, as illustrated by the model validation results for P. cattleianum at the 200 and 890 m sites. Heavily ¯owering trees may in¯uence the population time course of ¯owering, particularly if their ¯owering thresholds are similar, or if they are early or late ¯owering. To estimate these ¯owering thresholds and to ®x them in the simulation process is one way to overcome this problem. The general model [eqn (3)] and its different levels of simpli®cation can therefore respond to the different objectives and needs of a study, with different levels of simulation accuracy in consequence. The cost of higher accuracy is more precise observations, more complex parameter estimation and simulation limited to the ¯owering upon which parameters are ®tted. The ¯owering model can be connected to classical budburst models which provide the development time scale (e.g. Chuine et al., 1999). The critical development threshold used in the budburst models corresponds to the shoot or the tree ¯owering threshold of our model where it then follows a speci®ed distribution whose parameters are ®tted from the data. The model calculates the probability that a ¯ower bud bursts between ti ± 1 and ti. If individual ¯owers last 1 d, the model simulates directly the population of open ¯owers each day. But if individual ¯owers last nf days, the model simulates the onset of ¯owering for each ¯ower. The population of ¯owers open at ti is then calculated by adding the number of ¯owers that opened on the nf ± 1 days preceding ti (and which are still open) to the number of ¯owers that open between ti ± 1 and ti. Similar calculations can be made for particular events during the ¯owering

period, such as the period of pollen release or stigma receptivity (e.g. O'Brien and Calder, 1993), or the period in which ¯owers are susceptible to pests or pathogens. The ¯owering model has been built for and tested on tree species. However, it is applicable to any ¯owering species, including annual and herbaceous species. Requisites are the identi®cation of ¯owering units, as de®ned in the model construction, and the determination of a development time scale. Application of the model to Psidium cattleianum

Flowering in P. cattleianum shows both between- and within-tree asynchrony, the latter being more pronounced than the former. The `shoot' model is then the most relevant to describe the time course of ¯owering, as indicated by its high likelihood, whereas the `population' model has the lowest likelihood. Therefore, taking asynchrony into account in a ¯owering model improves the accuracy of the time course of ¯owering ®t for species affected by this phenomenon. Data show that ¯owering is not a smooth process at the tree or the population level (Figs 2±4). The irregularities are due to the simultaneous ¯owering of several shoots. The aggregative structure of the `shoot' and the `tree' models, and the integration of ¯owering variability, allow simulation of such irregularities (Fig. 1B and C), whereas the population model does not as it considers the ¯ower bud population as a whole and therefore smoothes the ¯owering process. The model parameters d and bs estimated using the 520 m data set, embrace the genetic and environmental variation of ¯owering phenology. Although further studies are needed to con®rm this result, the assumption of a similar normal distribution of the differences dj,p on all the trees is acceptable, indicating that the shoot ¯owering threshold distribution around the tree ¯owering threshold is not affected by the tree, and also that the different ¯owering shoots of a tree are comparable with respect to ¯owering. The tree ¯owering thresholds follow a normal distribution at the 520 m site. This distribution also ®ts the distribution of the tree ¯owering thresholds estimated at the sites at 200 and 890 m (n = 8, ks = 0´3045, P = 0´37; and n = 10, ks = 0´1873, P = 0´81, respectively). The genetic and environmental variability in tree ¯owering thresholds included in the normal distribution estimated at 520 m therefore appears representative of this variability at other sites with different genetic and environmental conditions. The hypothesis of a unique parameter of variance relative to ¯owering at the shoot level is not veri®ed. This parameter is highly variable between the ¯owering shoots of a tree, suggesting the effect of physiological or micro-environmental factors. Moreover, the distributions of this parameter value are signi®cantly different among individuals, indicating a tree effect, probably related to genetic and/or environmental factors. However, our data cannot be used to test this. Further work is needed to specify the factors affecting the parameter variability within and between trees. From a practical point of view, we retain this hypothesis for strawberry guava as the `shoot' model gives a satisfactory ®t of the data.

Normand et al. Ð Modelling Asynchronous Flowering Model applications

The model applications are wide and cover ®elds of research where simulation of the time course of ¯owering is useful: ¯owering phenology (e.g. Primack, 1980; Augspurger, 1983; Carthew, 1993); interactions between ¯oral biology and reproductive ecology [effective mating population; gene exchange (e.g. Fripp et al., 1987); outcrossing rates (e.g. Murawski and Hamrick, 1991, 1992); foraging behaviour of pollinators (e.g. Bawa, 1983); seed and fruit set (e.g. Carthew, 1993)]; distribution of airborne pollen release; integrated management of ¯ower pests and disease; and integration in a production model for crops whose pollination is a key step in seed and fruit set (Lescourret et al., 1999). The model is also a tool to quantify, at the individual or the population level, the within-tree variability and the tree ¯owering thresholds. These parameters are useful as selection criteria (Hof et al., 1999), to study variability in a population, or to estimate the contribution of genetic and environmental factors to ¯owering variability, given a relevant experimental design. L I TE R A T U R E C I T E D Agostini D, Habib R, Chadoeuf J. 1999. A stochastic approach for a model of ¯owering in kiwifruit `Hayward'. Journal of Horticultural Science and Biotechnology 74: 30±38. Andersen TB. 1991. A model to predict the beginning of the pollen season. Grana 30: 269±275. Augspurger CK. 1983. Phenology, ¯owering synchrony, and fruit set of six neotropical shrubs. Biotropica 15: 257±267. Bawa KS. 1983. Patterns of ¯owering in tropical plants. In: Jones CE, Little RJ, eds. Handbook of experimental pollination. New York: Scienti®c and Academic Editions, 394±410. Belmonte J, Roure JM. 1991. Characteristics of the aeropollen dynamics at several localities in Spain. Grana 30: 364±372. Boudry P, Wieber R, Saumitou-Laprade P, Pillen K, Van Dijk H, Jung C. 1994. Identi®cation of RLFP markers closely linked to the bolting gene B and their signi®cance for the study of the annual habit in beets (Beta vulgaris L.). Theoretical and Applied Genetics 88: 852±858. Boyer WD. 1973. Air temperature, heat sums, and pollen shedding phenology of longleaf pine. Ecology 54: 420±426. Cannell MGR, Smith RI. 1983. Thermal time, chill days and prediction of budburst in Picea sitchensis. Journal of Applied Ecology 20: 951± 963. Carthew SM. 1993. Patterns of ¯owering and fruit production in a natural population of Banksia spinulosa. Australian Journal of Botany 41: 465±480. Chuine I, Cour P, Rousseau DD. 1998. Fitting models predicting dates of ¯owering of temperate-zone trees using simulated annealing. Plant, Cell and Environment 21: 455±466. Chuine I, Cour P, Rousseau DD. 1999. Selecting models to predict the timing of ¯owering of temperate trees: implications for tree phenology modelling. Plant, Cell and Environment 22: 1±13. Citadin I, Raseira MCB, Herter FG, Baptista da Silva J. 2001. Heat requirement for blooming and lea®ng in peach. HortScience 36: 305± 307. DacuÈnha-Castelle D, Du¯o M. 1982. ProbabiliteÂs statistiques. I. ProbleÁmes aÁ temps ®xe. Paris: Masson. Dennis B, Kemp KP, Beckwith RC. 1986. Stochastic model of insect phenology: estimation and testing. Environmental Entomology 15: 540±546. Dodd JC, Estrada A, Jeger MJ. 1992. Epidemiology of Colletotrichum gloeosporioides in the tropics. In: Bailey JA, Jeger MJ, eds. Colletotrichum biology, pathology and control. Wallingford: CAB International.

415

Frenguelli G, Bricchi E, Romano B, Mincigrucci G, Spieksma FTM. 1989. A predictive study on the beginning of the pollen season for Gramineae and Olea europea L. Aerobiologia 5: 64±70. Fripp YJ, Grif®n AR, Moran GF. 1987. Variation in allele frequencies in the outcross pollen pool of Eucalyptus regnans F. Muell. throughout a ¯owering season. Heredity 59: 161±171. Gentry AH. 1974. Flowering phenology and diversity in tropical Bignoniaceae. Biotropica 6: 64±68. HaÈnninen H. 1987. Effects of temperature on dormancy release in woody plants: implications of prevailing models. Silva Fennica 21: 279±299. Hof L, Keiser LCP, Elberse IAM, Dolstra O. 1999. A model describing the ¯owering of single plants, and the heritability of ¯owering traits of Dimorphotheca pluvialis. Euphytica 110: 35±44. Ickovic MR, Boussioud-CorbieÁres F, Sutra JP, Thibaudon M. 1989. Hay fever symptoms compared to atmospheric pollen counts and ¯oral phenology within Paris suburban area in 1987 and 1988. Aerobiologia 5: 30±36. Ims RA. 1990. The ecology and evolution of reproductive synchrony. Tree 5: 135±140. Janick J, Moore JN. 1975. Advances in fruit breeding. West Lafayette: Purdue University Press. Kramer K. 1994. Selecting a model to predict the onset of growth of Fagus sylvatica. Journal of Applied Ecology 31: 172±181. Lehman EL. 1983. Theory of point estimation. New York: Wiley and Sons. Lescourret F, Blecher N, Habib R, Chadoeuf J, Agostini D, Pailly O, VaissieÁre B, Poggi I. 1999. Development of a simulation model for studying kiwi fruit orchard management. Agricultural Systems 59: 215±239. Mathsoft. 1999. S-Plus 2000 guide to statistics, volume 1. Seattle: Mathsoft. Medan D, Bartoloni N. 1998. Fecundity effects of dichogamy in an asynchronically ¯owering population: a genetic model. Annals of Botany 81: 373±383. Mitchell-Olds T. 1996. Genetic constraints on life-history evolution: quantitative-trait loci in¯uencing growth and ¯owering in Arabidopsis thaliana. Evolution 50: 140±145. Mosseler A, Papadopol CS. 1989. Seasonal isolation as a reproductive barrier among sympatric Salix species. Canadian Journal of Botany 67: 2563±2570. Murawski DA, Hamrick JL. 1991. The effect of the density of ¯owering individuals on the mating systems of nine tropical tree species. Heredity 67: 167±174. Murawski DA, Hamrick JL. 1992. Mating system and phenology of Ceiba pentandra (Bombacaceae) in Central Panama. Heredity 83: 401±404. Murfet IC. 1977. Environmental interaction and the genetics of ¯owering. Annual Review of Plant Physiology 28: 253±278. Murray MB, Cannell MGR, Smith RI. 1989. Date of budburst of ®fteen tree species in Britain following climatic warming. Journal of Applied Ecology 26: 693±700. Nienstaedt H. 1974. Genetic variation in some phenological characteristics of forest trees. In: Lieth H, ed. Phenology and seasonality modelling. New York: Springer Verlag, 389±400. Normand F, Habib R. 2001. Phenology of strawberry guava (Psidium cattleianum) in ReÂunion Island. Journal of Horticultural Science and Biotechnology 76: 541±545. O'Brien SP, Calder DM. 1993. Reproductive biology and ¯oral phenologies of the sympatric species Leptospermum myrsinoides and L. continentale (Myrtaceae). Australian Journal of Botany 41: 527±539. Osawa A, Shoemaker CA, Stedinger JR. 1983. A stochastic model of balsam ®r bud phenology utilizing maximum likelihood parameter estimation. Forest Science 29: 478±490. Pors B, Werner PA. 1989. Individual ¯owering time in a goldenrod (Solidago canadensis): ®eld experiment shows genotype more important than environment. American Journal of Botany 76: 1681±1688. Primack RB. 1980. Variation in the phenology of natural populations of montane shrubs in New Zealand. Journal of Ecology 68: 849±862. Sprent P. 1989. Applied nonparametric statistical methods. London: Chapman & Hall.