Andrea Veglio

The shade avoidance syndrome: A non-markovian stochastic growth model

Plants at high population density compete for light, showing a series of physiological responses known as the shade avoidance syndrome. These responses are controlled by the synthesis of the hormone auxin, which is regulated by two signals, an environmental one and an internal one. Considering that the auxin signal induces plant growth after a time lag, this work shows that plant growth can be modelled in terms of an energy-like function extremization, provided that the Markov property is not applied. The simulated height distributions are bimodal and right skewed, as in real community of plants. In the case of isolated plants, theoretical growth dynamics and speed correctly fit Arabidopsis thaliana experimental data reported in literature. Moreover, the growth dynamics of this model is shown to be consistent with the biomass production function of an independent model. These results suggest that memory effects play a nonnegligible role in plant growth processes.

Introduction

Ecologists have been concerned with competition in plants since the beginning of the last century. [START_REF] Clements | Research methods in ecology[END_REF] published the first formal definition of competition in plants, defining it as "the relation between plants occupying the same area and dependent upon the same supply of physical factors", and later [START_REF] Clements | Plant competition: an analysis of community functions[END_REF] he considered that competition begins "when the immediate supply of a single factor necessary falls below the combined demand of the individual plant".

Since the 50s a number of authors have systematically studied the effect of competition in plant systems. For example, [START_REF] Donald | Competition among pasture plants. i. intraspecific competition among annual pasture plants[END_REF] investigated in a seminal work the intra-specific light competition in pastures plants. More recently, [START_REF] Pacala | Forest models defined by field measurements: Estimation, error analysis and dynamics[END_REF] report experimental data providing little evidence of competition for nitrogen and water, and suggesting that the most important factor determining competition in sufficiently manured and watered plants is light. Accordingly, the present work will assume no stress due to insufficiency of nutrients or water (for models on dryland vegetation, see [START_REF] Kletter | Patterned vegetation and rainfall intermittency[END_REF] and ref.s therein) and will study light competition in plants grown at high population density, i.e. experiencing density stress. In particular, it will be concerned with intraspecific (inter-plant) light competition. The series of physiological responses demonstrated by plants experiencing density stress is collectively known as the shade avoidance syndrome and comprises effects such as stem elongation at Email addresses: andrea.veglio@unito.it (Andrea Veglio) expense of leaf and storage organ expansion, inhibition of branching, and acceleration of flowering [START_REF] Tao | Rapid synthesis of auxin via a new tryptophan-dependent pathway is required for shade avoidance in plants[END_REF].

There are two main 'fingerprints' characterizing the height distribution of a community of plants under density stress:

i. bimodality: bigger individuals grow more than smaller ones, thus become bigger and bigger, whilst smaller individuals grow much less [START_REF] Ford | Competition for light in a plant monoculture modelled as a spatial stochastic process[END_REF].

ii. hierarchy of exploitation (right skew): there is a large number of plants slightly smaller than the mean, and a small number of plants much greater than the mean, i.e. the distribution is right (positively) skewed [START_REF] Harper | A darwinian approach to plant ecology[END_REF].

The present work is based on the consideration that growth responses in plants at high population density are controlled by the synthesis of several phytohormones, the most important being auxin [START_REF] Friml | Plant biology: in their neighbour's shadow[END_REF][START_REF] Teale | Auxin in action: signalling, transport and the control of plant growth and development[END_REF]. Auxin induces plant growth after some delay (Taiz and Zeiger, 2002), and recent experimental results [START_REF] Tao | Rapid synthesis of auxin via a new tryptophan-dependent pathway is required for shade avoidance in plants[END_REF][START_REF] Stepanova | Taa1-mediated auxin biosynthesis is essential for hormone crosstalk and plant development[END_REF] show that the auxin pathway is mainly modulated by two signals, an environmental one and an internal one (or, using the vocabulary of Guédon et al. (2007), ontogenetic). The environmental signal consists in the decrease in the red to far-red ratio of incoming light. In fact, plants grown at high population density perceive both absorption of red light by canopy leaves and reflection of far-red light from neighbouring plants [START_REF] Tao | Rapid synthesis of auxin via a new tryptophan-dependent pathway is required for shade avoidance in plants[END_REF]. In other words, a plant 'recognizes' to be surrounded by other shading neighbours and responds producing auxin which induces growth after some time lag. On the other hand, the internal signal is mediated by the hormone ethylene (Stepanova et al A countless number of models have been developed to study both inter-and intra-specific competition among plants. Reviewing the different spatial community dynamics models, [START_REF] Bolker | Spatial dynamics in model plant communities: what do we really know?[END_REF] group them in three different frameworks, namely:

a. Interacting particle systems (IPS), which are stochastic, markovian [START_REF] Neuhauser | Mathematical challenges in spatial ecology[END_REF], continuous time models with discrete individuals located in the cells of a regular lattice (see, for example, [START_REF] Diggle | A spatial stochastic model of inter-plant competition[END_REF]; [START_REF] Durrett | Stochastic spatial models: A user's guide to ecological applications[END_REF]; [START_REF] Hendry | A coupled map lattice model of the growth of plant monocultures[END_REF]).

b. Stochastic point processes, which track discrete individuals that compete locally, but assume that they occupy a single point in continuous space (see, for example, [START_REF] Gandhi | critical slowing down" in time-to-extinction: an example of critical phenomena in ecology[END_REF]; [START_REF] Bolker | Spatial moment equations for plant competition: Understanding spatial strategies and the advantages of short dispersal[END_REF]; [START_REF] Dieckmann | The Geometry of Ecological Interactions: Simplifying Spatial Complexity[END_REF]).

c. Patch models, which track group of individuals, allowing (noncontiguous patch models) or not (metapopulation or patch-occupancy models) multiple individuals and multiple species per patch (see, for example, [START_REF] Levins | Regional coexistence of species and competition between rare species[END_REF]; [START_REF] Gilad | Ecosystem engineers: from pattern formation to habitat creation[END_REF][START_REF] Gilad | A mathematical model of plants as ecosystem engineers[END_REF]).

Strengths and weaknesses of each of the above frameworks have been exhaustively analyzed by [START_REF] Bolker | Spatial dynamics in model plant communities: what do we really know?[END_REF].

In this paper, I propose a model that studies the old problem of the light competition in terms of the old method of the energy-like function extremization, but does not belong to any of the previous frameworks. In fact, even sharing some features with other IPS models, the model presented here is non-markovian.

Indeed, markovianity is a very nice property in stochastic systems. It assures that in a system the future state depends only on the present state and is independent of past states [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF]. In other words, the system has no memory and this considerably simplifies calculations. Unfortunately, auxin production induces plant growth after a time lag, and at the proper time scale this memory effect forbids modelling the system by exploiting the Markov property.

I will show that in the present model the nonmarkovianity is necessary:

• to reproduce plant height distributions that are bimodal;

• to correctly mimic experimental growth dynamics in isolated plants;

• to be consistent with the biomass production model by [START_REF] Cournède | Computing competition for light in the greenlab model of plant growth: a contribution to the study of the effects of density on resource acquisition and architectural development[END_REF].

I will consider the total amount of auxin present in a plant, but not its spatial distribution in plant's body, i.e. the well characterized phenomenon of the auxin gradient (see Tanaka et al. (2006) and ref.s therein). For recent models on the auxin gradient the reader is referred to [START_REF] Jönsson | An auxin-driven polarized transport model for phyllotaxis[END_REF][START_REF] De Reuille | Computer simulations reveal properties of the cell-cell signaling network at the shoot apex in arabidopsis[END_REF]; [START_REF] Smith | A plausible model of phyllotaxis[END_REF]; [START_REF] Newell | Phyllotaxis: cooperation and competition between mechanical and biochemical processes[END_REF].

The model

As mentioned in the introduction, auxin is the main regulator of plant growth. For plants growing at high population density, auxin synthesis is regulated by: 1) an environmental signal due to the absorption of red light by canopy leaves and reflection of far-red light from neighboring plants [START_REF] Tao | Rapid synthesis of auxin via a new tryptophan-dependent pathway is required for shade avoidance in plants[END_REF];

2) an internal signal, which makes plants growing also when they do not experience density stress [START_REF] Stepanova | Taa1-mediated auxin biosynthesis is essential for hormone crosstalk and plant development[END_REF].

Consider a community of plants placed on a 2dimensional N × N square lattice. The lattice has N 2 nodes, each corresponding to a plant labelled as (i, j), with i, j = 1, . . . , N. Each (i, j)-plant is associated to its own amount of auxin s ij (t) at time t. Each plant shades, i.e. interact with, its neighbours by which conversely it is shaded. I will consider only the interaction between firstnearest neighbouring plants. Indeed, trials with secondnearest neighbour interactions showed that the results do not change substantially.

The auxin amount of a plant is converted into biomass after a delay (Taiz and Zeiger, 2002) that I will call α: if a certain auxin amount is produced by a plant at time t, then at time t + α the plant shoots to grow proportionally to the produced quantity of auxin (in this paper I will take into account only vertical shooting, i.e. I will assume the height of each plant to be proportional to its weight.) Accordingly, I will define the shade avoidance potential a ij (t) as the function representing the predisposition at time t of the (i, j)-plant to avoid shade by shooting at time t + α. Its definition needs some steps.

For simplicity, first consider an ideal subsystem of two plants, say the (i, j)and the (i + 1, j)-plant, subjected only to the environmental signal 1).

The predisposition at time t of the (i, j)-plant to shoot at time t + α depends first of all on its own auxin amount s ij (t). In fact, the greater s ij (t) is, the more the (i, j)plant has the possibility to grow at time t + α.

Moreover, the (i, j)-plant perceives an environmental signal from the (i + 1, j)-plant, that absorbs red light and reflects far-red light. This signal "warns" the (i, j)-plant of the presence of its (i + 1, j) neighbour. The greater the (i + 1, j) auxin amount s i+1j (t) is at time t, the more (i + 1, j) will shade the (i, j)-plant at time t + α. Because the (i, j)-plant "wants" to avoid the (i + 1, j) shade, its predisposition at time t to shoot at time t + α has to depend also on the auxin amount s i+1j (t) of its neighbour. In fact, the greater s i+1j (t) is, the more the (i, j)-plant has to grow to avoid the shade of the (i + 1, j)-plant at time t + α.

Therefore, for a system of two plants subjected only to the environmental signal, the predisposition at time t of the (i, j)-plant to shoot at time t + α, i.e. the shade avoidance potential a i,j (t), reads

a ij (t) = J s ij (t) • s i+1j (t).
(1)

A c c e p t e d m a n u s c r i p t

J is a non-negative constant that expresses how strong plants shade each other or, in other words, how strong they have to compete for light. Because in a generic system of N2 elements the (i, j)-plant has more than one firstnearest neighbour (namely, two if (i, j) is at the corner of the lattice, three if it is at the border, four if inside), the (i, j)-plant shades and is shaded by all its first nearest neighbours. Therefore, all the contributions from the firstnearest neighbour have to be taken into account in the shade avoidance potential, i.e.

a ij (t) = J s ij (t) n.n. s ij (t), (2)
where

n.n. s ij (t) ≡ s i-1j (t)+s i+1j (t)+s ij-1 (t)+s ij+1 (t)
is the sum over the first nearest neighbour of the (i, j)plant.

Eq. ( 2) takes into account only the environmental signal 1) perceived by the (i, j)-plant; if it had no neighbours, its shade avoidance potential would be zero, i.e. it would have no chance to produce biomass. Since isolated plants grow, the internal signal 2) can be taken into account by adding a term proportional to the auxin amount in the plant, hs ij (t), where h is a non-negative constant expressing the strength of the signal, i.e. how high the plant would grow if isolated. Therefore, the complete function expressing the shade avoidance potential of the generic (i, j)-plant reads

a ij (t) = h s ij (t) + J s ij (t) n.n.
s ij (t).

(3)

Notice that the greater the shade avoidance potential is at time t, the more the (i, j)-plant will be able to shoot at time t + α, i.e. the more it will be able to intercept light and avoid shade.

Studying the system in the mean field approximation, it is worth considering the global shade avoidance potential A(t) ≡ N i=1 N j=1 a ij (t). Because all the plants of the community contribute with their mutual shading, A(t) is a global quantity that represents the predisposition of the whole system to avoid shade by shooting after the time lag or, in other words, the mean shade avoidance potential up to a N -2 factor. A(t) reads

A(t) = ij s ij (t) h + J n.n. s ij (t) , (4)
where the sum ij ≡ N i=1 N j=1 runs over all the N 2 plants. In passing, notice that Eq. ( 4) is formally analogous to the Ising hamiltonian1 [START_REF] Huang | Statistical Mechanics[END_REF].

Accordingly to what discussed above, it is reasonable to assume that on average plants try to maximize the shade avoidance potential: the greater is the global potential at time t, the maximal will be the light exposure of the community at time t + α. Therefore, at each time step t the auxin amount of a randomly selected (i, j)-plant may be increased only if the global shade avoidance potential is maximized, as follows:

s ij (t + 1) = s ij (t) ( 5 ) + ρ ξ(t) e -bt Θ A(t) -A(t -α) . ξ(t) ∈ [-1, 1]
is a randomly drawn number with flat distribution, ρ > 0 the maximum updating amplitude for the selected auxin amount, 1/b > 0 the equilibrium time 2 . Θ[•] is the Heaviside step function, defined as

Θ[x] = 0 if x ≤ 0 and Θ[x] = 1 if x > 0.
The exponential term e -bt has been chosen to let the growing process dramatically slow down once reached the equilibrium time 1/b. The positive discrete variable α represents the time lag introduced above and is calculated according to the flowchart in Fig. 1. Because at the very beginning of a plant life (dormant seed) the auxin amount is negligible, it is reasonable to choose as initial condition s ij (0) = 0.

Eq. ( 5) states that, at each time step, the selected auxin amount s ij can be updated by a randomly drawn quantity (which gets smaller and smaller as time goes by because of the exponential term) only if such an update implies the maximization (Θ[A(t) -A(tα)]) of the shade avoidance potential (4), whose dynamics is defined as:

A(t + 1) = A(t) ( 6 )
+ ρ ξ(t) e -bt h + J n.n.

s ij (t) ,

with A(0) = 0. ξ(t) is the same random number selected in (5). The sum n.n. runs over the nearest neighbours of the selected (i, j)-plant as above.

A mortality rate has been introduced in the following way: if at certain time the (i, j)-plant has reached an auxin amount smaller than a threshold d, s ij is set to zero and no longer updated during the simulation. This corresponds to the death of the (i, j)-plant. For the sake of clarity, the mortality process has not been shown in the flowchart in Fig. 1. Now, consider the function q ij (t), that shall represent the height of the (i, j)-plant at time t. As discussed before, the auxin amount of a plant is converted into biomass after a time lag α, and the biomass of a plant is assumed to be proportional to its height, q ij (t). Therefore, setting the proportionality constant equal to one, before the equilibrium is reached the plant height at time t corresponds to the auxin amount at time tα. Then, at the equilibrium, the system does not evolve appreciably and, therefore, the plant height at time t corresponds to the auxin amount at the same time. In formulae,

q ij (t) = s ij (t -α) fort < 1/b, (7) q ij (t) = s ij (t) f o rt 1/b. ( 8

Numerical results

I performed numerical simulations of the dynamics (5-6) using the following parameter values: N = 10, ρ = .2, b = 10 -4 , J = 10 2 , h = 1. The mortality threshold has been set to d = .3, which corresponds to a rate of mortality of 70%, accordingly to experiments by [START_REF] Ford | Competition and stand structure in some evenaged plant monocultures[END_REF].

Fig. 2 shows two height distributions at the equilibrium, i.e. at t 1/b, for initial communities of N 2 = 100 plants. Each histogram represents the height distribution of the dN 2 = 30 survived plants.

In Fig. 2a the simulation has been performed according to the model presented in the previous section: the system is non-markovian. The distribution results to be bimodal and right-skewed (see experiments by [START_REF] Ford | Competition and stand structure in some evenaged plant monocultures[END_REF]). Moreover, plants at the lattice border are on average (see Appendix B) shorter than central ones.

On the contrary, in Fig. 2b, the time lag α has been set to be equal to zero, i.e. the system has been forced to be markovian. The height distribution looses bimodality, showing that in this model framework the nonmarkovianity is a fundamental feature to interpret data of plants competing for light.

The numerical simulations show that α varies with time. In particular, α can be consistently fitted by the simple saturating function

α(t) = m t K + t , ( 9
)
where m is the plateau value reached by α and K ≡ 1/(2b).

Isolated plants

The model presented here is intended to mimic the behaviour of plants at high population density. Indeed, it may be extended to isolated plants. 1/b, for an initial community of N 2 = 100 plants. The rate of mortality have been set to 70%, accordingly to experiments by [START_REF] Ford | Competition and stand structure in some evenaged plant monocultures[END_REF]. The heights are reported in percentage relatively to the maximal height in each sample. a) Simulation of the system according to the model: the Markov property is relaxed. The distribution is bimodal and right-skewed, as in the experiments by [START_REF] Ford | Competition and stand structure in some evenaged plant monocultures[END_REF]. b) Simulation of the system setting the time lag α = 0, i.e. forcing the system to be Markovian. The distribution looses bimodality.

A c c e p t e d m a n u s c r i p t

In Appendix A, I derived the expression for the growth (height) dynamics q(t) and the growth speed q(t) in isolated plants. They read

q(t) = ρ Ω b -1 (1 -e -bf (t) ), (10)
q(t) = ρ Ω 1 - m K (K + t) 2 e -bf (t) , ( 11
)
where Ω is defined as in (A3) and

f (t) = t -α(t) = K -m + t K/t + 1 (12)
This allows to directly check the model consistency with data available in literature. [START_REF] Jouve | Arabidopsis thaliana floral stem elongation: Evidence for an endogenous circadian rhythm[END_REF] report experimental data for Arabidopsis thaliana floral stem elongation over time. They study circadian rhythms and report experimental data for Arabidopsis first inflorescence internode length and for first inflorescence extension rate as functions of time. The first inflorescence internode length can be assumed to be proportional to the plant growth dynamics q(t), and the first inflorescence extension rate, up to a constant, to the growth speed q(t).

Fig. 3 shows q(t) (10) fitted on Arabidopsis first inflorescence internode length data and Fig. 4 shows q(t) (11) fitted on Arabidopsis first inflorescence extension rate data.

The agreement are satisfactory in both cases; fluctuations in first inflorescence extension rate data are due to the circadian rhythm, which is not taken into account by this model and which can be easily averaged out.

Interestingly, [START_REF] Cournède | Computing competition for light in the greenlab model of plant growth: a contribution to the study of the effects of density on resource acquisition and architectural development[END_REF] propose the following empirical formula for the biomass production Q(t),

Q(t) = β λ(t) (1 -e -γS(t) ), ( 13
)
where β and γ are characteristic parameters, S(t) the leaf surface area of the plants, and λ(t) a function of environmental conditions related to the evapotranspiration that I here assume to be constant. Considering again the Arabidopsis case, experimental data from Fig. 3a in [START_REF] Cookson | Day length affects the dynamics of leaf expansion and cellular development in arabidopsis thaliana partially through floral transition timing[END_REF] show that, in the time window of interest (0-8 days), Arabidopsis leaf surface area over time can be consistently fitted (reduced χ 2 ∼ 1) by af (t), being a a proportionality constant.

This suggests that the growth dynamics q(t) obtained here and the biomass production Q(t) proposed by [START_REF] Cournède | Computing competition for light in the greenlab model of plant growth: a contribution to the study of the effects of density on resource acquisition and architectural development[END_REF] are consistent. While parameter identification between the two models needs further investigation, it is worth observing that f (t) can correctly reproduce experimental data for leaf surface area over time only in the case m = 0. In other terms, it is necessary the system to keep memory of the time lag α even in the case of isolated plants.

Sensitivity analysis

Focusing back to the general model (J ≥ 0, N ≥ 1), it is interesting to analyze its parameter sensitivity.

A c c e p t e d m a n u s c r i p t

I have performed a stochastic analysis of the model (Appendix B), whose main result is, according to (B13) and (B14),

s ij st function(ρ, b, N ). ( 14
)
This means that on average the auxin amount in a generic (i, j)-plant does not directly depend on J and h. Numerical simulations confirmed this result. Speculating, this suggests that plant height does not directly depend either on the R:FR ratio or on the hormone ethylene density. This observation is consistent with the fact that plant sensitivity both to light [START_REF] Maloof | Natural variation in light sensitivity of arabidopsis[END_REF] and to hormone ethylene [START_REF] Chen | Analysis of ethylene signaltransduction kinetics associated with seedling-growth response and chitinase induction in wild-type and mutant arabidopsis[END_REF] stimulations rapidly saturates.

Discussion

This paper proposes a simple model for the growth at high population density of plants competing for light. The model is based on the consideration that the plant growth is controlled by the auxin synthetic pathway, which is regulated by two signals, an internal one related to the hormone ethylene, and an environmental one, due to the mutual shading. Moreover, the auxin signal is converted into biomass only after a time lag, α, that varies on time. This memory effect forbids modelling the system markovianly at the natural time scale. In fact, a single plant has the predisposition to avoid the shade experienced at time t only by shooting at time t + α. At the system level, this can be modelled by maximizing the global shade avoidance potential (4).

Consequently, this model cannot be considered to belong to any of the three frameworks of spatial community models drafted by [START_REF] Bolker | Spatial dynamics in model plant communities: what do we really know?[END_REF], namely the Interacting Particle Systems, the Stochastic Point Processes, and the Patch Models. In fact, even sharing some features with IPSs, the model proposed here is non-markovian. Indeed, it can be considered closer to those models that describe in terms of energy-like function(al)s a broad range of life science phenomena, such as X chromosome inactivation [START_REF] Nicodemi | Self-assembly and dna binding of the blocking factor in x chromosome inactivation[END_REF], epithelial cell polarization [START_REF] Veglio | Symmetry breaking mechanism for epithelial cell polarization[END_REF], and eukaryotic chemotaxis [START_REF] Ferraro | Spatial signal amplification in cell biology: a lattice-gas model for selftuned phase ordering[END_REF][START_REF] Gamba | Universal features of cell polarization processes[END_REF].

The present work shows that the simulated plant height distributions are bimodal and right skewed, as in real communities of plants at high population density. On the contrary, if no delay is assumed between the auxin signal and the biomass production, i.e. if the time lag α is set to zero, simulations show that the distribution of plant heights is unimodal and cannot reproduce real distributions.

Extending the model to isolated plants, an analytical expression is obtained both for the growth dynamics (10) and the growth speed (11). They have been fitted to Arabidopsis thaliana first inflorescence length data and on first inflorescence internode extension rate data, respectively. Theoretical functions and experimental data show a good agreement.

Moreover, the obtained growth dynamics (10) is consistent with the biomass production function proposed by [START_REF] Cournède | Computing competition for light in the greenlab model of plant growth: a contribution to the study of the effects of density on resource acquisition and architectural development[END_REF].

The bimodality in height distribution, the good match between experimental data and theoretical functions for isolated plants, and the consistency with an empirical model independently derived, can be considered the main results of this work. They hold as long as the time lag α is taken into account. If the model is forced to be markovian, α = 0, real distributions and dynamics can no longer be mimicked. This suggests that memory effects play a key role in the plant growth process and that they cannot be neglected.

Interestingly, [START_REF] Hara | A stochastic model and the moment dynamics of the growth and size distribution in plant populations[END_REF] studies the dynamics of plant populations with a stochastic model based on the Kolmogorov forward equation, which needs markovianity, and observes that the Markov property applied to plant growth is reasonable as a first approximation and needs further investigation. In fact, it seems not to be a case that his model does not reproduce bimodality in plant height distributions.

A further step in the comprehension of plant growth would be to experimentally check the effects of varying the time lag between the auxin production and its conversion in biomass (for example, along the lines of [START_REF] Tao | Rapid synthesis of auxin via a new tryptophan-dependent pathway is required for shade avoidance in plants[END_REF]), both in high population density and isolated plants. Accordingly, from a theoretical point of view it would be interesting to study the bimodal-to-unimodal switch in height distribution due to the loss of markovianity, in particular in terms of statistical physics of phase transitions [START_REF] Huang | Statistical Mechanics[END_REF][START_REF] Barabási | Fractal Concepts in Surface Growth[END_REF].

A c c e p t e d m a n u s c r i p t

which has the following steady state,

s ij st ≡ lim t→∞ s ij (t) = ρ Ω ij b N 2 . (A6)
The dynamics of an isolated plant height over time q(t) may be defined as

q(t) ≡ s(t) ⏐ ⏐ t=t-α(t) , (A7)
where ṡ(t) is the same as in (A4), but without position indexes because it is referred to an isolated plant. It is calculated at time tα(t) to reflect that auxin is not immediately converted into biomass, as discussed in Section 2. (A7) is the same as ( 7) except for the average, which is intended to reduce fluctuations due to the isolated plant condition. The shade avoidance potential is ideally evaluated in the condition J = 0, i.e. no interaction, and N = 1. Using ( 9), ( A5) and (A7), q(t) reads

q(t) = s(t) ⏐ ⏐ t=t-α(t) (A8) = ρ Ω b -1 (1 -e -b(t-mt/(K+t)) ), and
q(t) = ρ Ω 1 - m K (K + t) 2 e -b(t-mt/(K+t)) . (A9)
The fact that α tends to the finite value m implies that at the equilibrium the delay gets negligible, letting q(t) = s(t) (A5) and q(t) = ṡ(t) (A4). This is consistent with (8), that holds after the equilibrium.

B. Analysis of the Parameters

In order to explain the model dependence on the parameters, it is worth evaluating analytically the noise mean Ω ij introduced in (A3). Its definition is general and valid for the whole system of plants.

(A1) is a Langevin-like equation with a non-zero-mean colored noise, no drift term and time dependent diffusion term. It is possible to perform some transformations to obtain a more standard Langevin-like equation. Defining

x ij (t) ≡ s ij (t)e bN 2 t , (B1)
and inserting (A1) in the time derivative of (B1), I obtain

ẋij (t) = F ij (x) + ρ ηij (t), (B2)
where

F ij (x) ≡ bN 2 x ij (t) + ρΩ ij , (B3) ηij (t) ≡ η ij (t) -Ω ij . (B4)
Consider that x(t) and s(t) have, at the steady state, the same pdf.

(B2) is a Langevin-like equation with time independent drift term, constant diffusion term and zero-mean colored noise.

Assuming that such a noise is a Ornstein-Uhlenbeck process (van Kampen, 2007), i.e.

dη ij (t) dt = - 1 μ ij ηij (t) + ζ(t), (B5)
where μ ij is the correlation time and ζ(t) is a white noise with zero mean and delta time-correlation, it is possible to write the Fokker-Planck equation corresponding to (B2), i.e.

∂p ij (x, η, t) ∂t = -bN 2 + γp ij (x, η, t) + (B6) + γη ∂p ij (x, η, t) ∂η + 1 2 ∂ 2 p ij (x, η, t) ∂η 2 .
Even at the steady state, (B6) is not easy to be analytically solved. However, for a huge set of parameters (B2) lies in the validity range of the Unified Colored Noise Approximation, UCNA [START_REF] Jung | Dynamical systems: A unified colorednoise approximation[END_REF][START_REF] Hänggi | Ch. Colored Noise in Dynamical Systems[END_REF][START_REF] Hasegawa | Dynamics of the langevin model subjected to colored noise: Functional-integral method[END_REF]. In fact, to apply the UCNA it is necessary [START_REF] Hänggi | Ch. Colored Noise in Dynamical Systems[END_REF] that

x ij μ ij b . (B7)
According to (B1) and to the fact that s ij ∼ 1, constraint (B7) corresponds to

t log(μ ij /b) bN 2 ≡ t u . (B8)
Therefore, since equilibrium time 1/b t u , the UCNA can be applied for the parameter set used in the reported numerical simulations. Applying UCNA, the following effective Langevin equation holds

ẋij (t) σ ij F ij (x) + ρ σ 2 ij ζ(t), (B9)
where ζ(t) is again a white noise with zero mean and delta time-correlation, while I defined σ ij as

σ ij ≡ 1 1 -μ ij bN 2 .
(B10) I calculated the correlation time μ ij by averaging the numerical values of α over all the lattice nodes and over different realizations of the system. As an example, labelling the site at the center of the lattice as cc, the sites at the center of the lattice borders as bb and the sites at the corners as aa, I obtained, in continuous time units, μ cc 54, μ bb 50 and μ aa 47.

In passing, notice that the drift term of (B9), σ ij F ij (x), is linear in x ij ; therefore this model lies in a class of universality different from the principal growth models3 [START_REF] Barabási | Fractal Concepts in Surface Growth[END_REF]

  which stimulates the production of auxin even in isolated plants.

Figure 1 :

 1 Figure 1: Flowchart of the system dynamics.

Figure 2 :

 2 Figure2: Height distribution at the equilibrium, t 1/b, for an initial community of N 2 = 100 plants. The rate of mortality have been set to 70%, accordingly to experiments by[START_REF] Ford | Competition and stand structure in some evenaged plant monocultures[END_REF]. The heights are reported in percentage relatively to the maximal height in each sample. a) Simulation of the system according to the model: the Markov property is relaxed. The distribution is bimodal and right-skewed, as in the experiments by[START_REF] Ford | Competition and stand structure in some evenaged plant monocultures[END_REF]. b) Simulation of the system setting the time lag α = 0, i.e. forcing the system to be Markovian. The distribution looses bimodality.

Figure 3 :

 3 Figure 3: Arabidopsis thaliana growth dynamics. Eq. (11) (continuous line) and experimental data from Fig. 1, dataset A (12 hours of dark every 12 hours of light), in Jouve et al. (1998) (dashed line).

Figure 4 :

 4 Figure 4: Arabidopsis thaliana growth speed. Eq. (10) (continous line) and experimental data from Fig. 2A in Jouve et al. (1998) (dashed line).

Figure 5 :

 5 Figure5: Eq. (A5) (lines) and numerical results from the simulation of the dynamics (5-6) (symbols). Arbitrary units.

  . Then, it is possible to write the Fokker-Planck equation corresponding to (B9),

	∂p ij (s, t) ∂t	= -σ ij	∂ F ij (s)p ij (s, t) ∂s	+
		+	1 2	ρ 2 σ 4 ij	∂ 2 p ij (s, t) ∂s 2	,	(B11)

Actually, it is the opposite. In fact, here the global shade avoidance potential will be maximized, while in the Ising model the hamiltonian is minimized.

In this work, for equilibrium I mean the condition when the increments Δs ij become negligible. This happens at time t ∼ 1/b.

In particular, in the Kardar-Parisi-Zhang model[START_REF] Kardar | Dynamic scaling of growing interfaces[END_REF] the drift term is equal to ν∇ 2 x ij + μ 2 (∇x ij ) 2 , where ν and μ are characteristic parameters of that model.

Acknowledgements

The author wishes to thank Federico Bussolino, Michele Caselle, Antonio Celani, and Marta Lucchini for useful discussions and comments, and Guido Serini for carefully reading the manuscript.

Appendices

A. Analysis of the Dynamics

To analyze the dynamics, it is worth switching from discrete to continuous time by substituting t → N 2 t. This yields (5) to read

where

is a random term. I will treat it as a gaussian non-zeromean colored noise with time-independent mean, which I shall label as

Therefore, the average of (A1) over different system's realizations is simply

with the straightforward solution )


Thibaut Jombart
email: tjombart@imperial.ac.uk

Sandrine Pavoine

Sébastien Devillard

Dominique Pontier

Putting phylogeny into the analysis of biological traits: A methodological approach

Keywords: phylogenetic principal component analysis, pPCA, autocorrelation, multivariate, comparative method, phylogenetic signal

Phylogenetic comparative methods have long considered phylogenetic signal as a source of statistical bias in the correlative analysis of biological traits. However, the main life-history strategies existing in a set of taxa are often combinations of life history traits that are inherently phylogenetically structured. In this paper, we present a method for identifying evolutionary strategies from large sets of biological traits, using phylogeny as a source of meaningful historical and ecological information. Our methodology extends a multivariate method developed for the analysis of spatial patterns, and relies on finding combinations of traits that are phylogenetically autocorrelated. Using extensive simulations, we show that our method efficiently uncovers phylogenetic structures with respect to various tree topologies, and remains powerful in cases where a large majority of traits are not phylogenetically structured. Our methodology is illustrated using empirical data, and implemented in the free software R.

A c c e p t e d m

a n u s c r i p t

Introduction

Phylogeny has long been recognised as a major source of biological variation.

For instance, [START_REF] Gregory | Convergence and applied phenomena in the mammalia[END_REF] and [START_REF] Osborn | Heritage and habitus[END_REF] considered that species' variability should be partitioned between heritage (i.e., phylogenetic inertia) and habitus (i.e., adaptation). In their well-known criticism of the adaptationist paradigm, [START_REF] Gould | The spandrels of san marco and the panglossian paradigm: a critique of the adaptationist program[END_REF] underlined the importance of the constraints imposed by the phylogeny to the variability observed among organisms.

In comparative studies, the effect of phylogeny has merely been perceived as a source of nuisance, since it reveals non-independence among trait values observed in taxa [START_REF] Dobson | The use of phylogeny in behavior and ecology[END_REF][START_REF] Faraway | Phylogenies and the comparative method[END_REF], and thus violates one of the basic assumptions required by most statistical tools [START_REF] Harvey | The Comparative Method in Evolutionary Biology[END_REF].

Phylogenetic comparative methods (PCM) were especially designed to solve this problem. Various methods have been developed that transform quantitative traits into new variables that are not correlated to phylogeny, according to a given model of evolution. For instance, phylogenetic independent contrasts (PIC, Felsenstein, 1985) transform values observed at the n tips of a phylogeny into (n -1) node values that are not phylogenetically autocorrelated under a Brownian motion model. Generalised least squares (GLS, [START_REF] Grafen | The phylogenetic regression[END_REF][START_REF] Rohlf | Comparative methods for the analysis of continuous variables: geometric interpretations[END_REF]) is a more general technique that allows specifying the autocorrelation of observations as a component of a linear model. This approach can therefore account for the non-independence among observations using a wide variety of models of evolution [START_REF] Hansen | Translating between microevolutionary process and macroevolutionary patterns: the correlation structure of interspecific data[END_REF]. As stressed by [START_REF] Rohlf | A comment on phylogenetic correction[END_REF], these approaches do not actually remove phylogenetic autocorrelation from the data, but merely take it into account to provide more accurate estimates of model parameters. In fact, PIC, GLS, along with other existing PCM all aim towards the same goal: 'correcting for phylogeny' in the correlative analysis of biological traits at the species level [START_REF] Harvey | Comparative methods for explaining adaptations[END_REF][START_REF] Martins | Adaptation and the comparative method[END_REF][START_REF] Martins | Adaptive constraints and the phylogenetic comparative method: a computer simulation test[END_REF][START_REF] Garland | Phylogenetic approaches in comparative physiology[END_REF].

A c c e p t e d m a n u s c r i p t

Nonetheless, studying the phylogenetic patterns of trait variation allows formation of hypotheses about the evolutionary pathways that led to the trait values of extant species. It also allows shedding light onto the influence of historical and ecological processes on community assembly [START_REF] Webb | Phylogenies and community ecology[END_REF].

Many biologically meaningful patterns are inherently structured with phylogeny.

Indeed, many life-history and ecological strategies are likely to be phylogenetically structured [START_REF] Webb | Phylogenies and community ecology[END_REF]. Inheritance from a common ancestor and phylogenetic inertia (i.e., constraints to evolution) may cause phylogenetic signal (similar trait values across closely related species) to occur. Other factors leading to phylogenetic signals in traits act at the population level rather than at the species level such as high gene flow, lack of genetic variation, stabilising selection if changes in trait states reduce fitness, or population growth if traits are pleiotropically linked to other traits that reduce fitness [START_REF] Wiens | Niche conservatism: integrating evolution, ecology, and conservation biology[END_REF]. However, traits might also be affected by variations unrelated to the phylogeny, but relating to ecological conditions experienced by the species.

For instance, biotic interactions might drive character displacement and abiotic interactions might lead to trait convergence. From this perspective, phylogenetic signal becomes a source of precious biological information that can be used to identify historical as well as recent evolutionary strategies. Interestingly, a similar paradigm shift occurred in spatial ecology [START_REF] Legendre | Spatial autocorrelation: trouble or new paradigm?[END_REF] when it was pointed out that spatial patterns in species' distribution were not only sources of spurious correlations, but also indicators of critical ecological structures such as localised species assemblages and species-environment associations. This paradigm shift proved particularly fecund and still motivates innovative developments in statistical ecology (e.g., [START_REF] Dray | Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM)[END_REF][START_REF] Griffith | Spatial modeling in ecology: the flexibility of eigenfunction spatial analyses[END_REF].

In this paper, we present a method which uses phylogenetic information to uncover the main phylogenetic structures observable in multivariate data associated with a phylogeny. Our approach, phylogenetic principal component analysis (pPCA), extends a methodology developed in spatial ecology (Dray , 2008) and in spatial genetics [START_REF] Jombart | Revealing cryptic spatial patterns in genetic variability by a new multivariate method[END_REF] to the analysis of phylogenetic structures in biological features of taxa such as life-history traits.

We emphasise that phylogenetic structures can be measured and quantified in the same way as spatial structures, as they are both associated with the concept of autocorrelation. We then define different kinds of phylogenetic structures, and show how pPCA can be used to identify them. After evaluating the ability of pPCA to uncover phylogenetic patterns through extensive simulations, we illustrate our method using an empirical example. pPCA is implemented in the adephylo package [START_REF] Jombart | adephylo: exploratory analyses for the phylogenetic comparative method[END_REF] for the free software R (R Development Core Team, 2009).

Methods

Measuring phylogenetic autocorrelation

Phylogenetic autocorrelation is said to occur whenever the values taken by a set of taxa for a given biological trait are not independent of the phylogeny. Frequently, closely related taxa exhibit more similar traits than randomly-chosen taxa. Moran's I, an index originally used to measure spatial autocorrelation [START_REF] Moran | The interpretation of statistical maps[END_REF][START_REF] Moran | Notes on continuous stochastic phenomena[END_REF], has been proposed for measuring phylogenetic autocorrelation [START_REF] Gittleman | Adaptation: statistics and a null model for estimating phylogenetic effects[END_REF]. Adapting the former definition (Cliff and Ord, 1973, p13) to the phylogenetic context, I is defined as:

I W (x) = x T Wx n 1 var(x) (1)
where x is the centred vector of a trait observed on n taxa, var(x) is the usual variance of x, and W is a matrix of phylogenetic proximities among taxa

(W = [w ij ] with i, j = 1, . . . , n), whose diagonal terms are zero (w ii = 0),
and rows sum to one ( n j=1 w ij = 1). The null value, i.e. the expected value when no phylogenetic autocorrelation arises, is [START_REF] Cliff | Spatial autocorrelation[END_REF]. In its initial formulation [START_REF] Gittleman | Adaptation: statistics and a null model for estimating phylogenetic effects[END_REF], i.e. before row standardisation so that n j=1 w ij = 1, W contained binary weights. Before this standardisation, the entry at row i and column j was set to 1 if taxon i shared

I 0 = -1/(n -1)

A c c e p t e d m a n u s c r i p t

a common ancestor with taxon j at a given taxonomic level, and to 0 otherwise.

Hence, taxa were considered as either phylogenetically related or not. Moran's I then compared the trait value of a taxon to the mean trait value in related taxa to detect phylogenetic autocorrelation.

Such binary relationships are clearly not sufficient to model the possibly complex structure of proximities among taxa induced by the phylogeny. To achieve better resolution in these comparisons, we propose using as entries of W any measurement of phylogenetic proximity valued in R + verifying:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ w ij 0 ∀ i, j = 1, . . . , n w ii = 0 ∀ i = 1, . . . , n n j=1 w ij = 1 ∀ i = 1, . . . , n (2)
Then, Moran's I compares the value of a trait in one taxon (terms of x)

to a weighted mean of other taxa states (terms of Wx) in which phylogenetically closer taxa are given stronger weights. This extension gives the index considerable flexibility for quantifying phylogenetic autocorrelation, as phylogenetic proximities can be derived from any model of evolution (including or not branch lengths). For instance, one interesting possibility would be using the covariance matrix estimated in a GLS model [START_REF] Grafen | The phylogenetic regression[END_REF] to define phylogenetic proximities. This could be achieved by setting diagonal terms (variances) of the covariance matrix to zero, adding the smallest constant ensuring that all terms are positive, and row-standardizing the resulting matrix.

This formulation of Moran's I also relates the index to other PCM. For instance, the test proposed by [START_REF] Abouheif | A method for testing the assumption of phylogenetic independence in comparative data[END_REF], initially based on the many possible planar representations of a tree, turned out to be a Moran's I test using a particular measure of phylogenetic proximity for W [START_REF] Pavoine | Testing for phylogenetic signal in life history variable: Abouheif's test revisited[END_REF]).

Moran's I is also related to autoregressive models. In their simplest form,

A c c e p t e d m a n u s c r i p t

these models are written as (Cheverud and Dow, 1985;Cheverud et al., 1985):

x = ρWx + Zβ + e (3)
where ρ is the autocorrelation coefficient, Z is a matrix of explanatory variables, β is the vector of coefficients, and e is a vector of residuals. The matrix of phylogenetic relatedness W (Cheverud and Dow, 1985;Cheverud et al., 1985) is exactly the weight matrix of our definition of Moran's I (equation 1). The essential difference between the two approaches is that autoregressive models perform the regression of x onto Wx, while I computes the inner product between both vectors (numerator of equation 1) to measure phylogenetic autocorrelation.

Lastly, the weighting matrix W is also the core of another approach producing variables that model phylogenetic structures [START_REF] Peres-Neto | A unified strategy for estimating and controlling spatial, temporal and phylogenetic autocorrelation in ecological models[END_REF]. Like

Moran's I, this approach was initially developed in spatial statistics [START_REF] Griffith | Spatial autocorrelation and eigenfunctions of the geographic weights matrix accompanying geo-referenced data[END_REF], and consisted in finding eigenvectors of a doubly centered spatial weighting matrix [START_REF] Dray | Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM)[END_REF]. Applied to a matrix of phylogenetic proximity W, this method yields uncorrelated variables modeling different observable phylogenetic patterns, each related to a value of Moran's I. [START_REF] Peres-Neto | A unified strategy for estimating and controlling spatial, temporal and phylogenetic autocorrelation in ecological models[END_REF] performed the regression of a variable x onto these eigenvectors to partial-out the phylogenetic autocorrelation from x. Alternatively, we suggest using these eigenvectors to simulate what we further call 'global' and 'local' phylogenetic structures. of taxa are not independent, but tend to be more similar in closely related taxa (e.g., Figure 1A). Most common explanations for this phenomenon are inheritance from a common ancestor, or the conservation of ecological niches [START_REF] Harvey | The Comparative Method in Evolutionary Biology[END_REF].Traits whose evolution can be modeled by a Brownian or by an OU process with low stabilising constraint generally display global patterns [START_REF] Abouheif | A method for testing the assumption of phylogenetic independence in comparative data[END_REF][START_REF] Pavoine | Testing for phylogenetic signal in life history variable: Abouheif's test revisited[END_REF]. Such phenomenon typically results in close-to-the-root divergence in evolutionary strategies.

Global and local phylogenetic structures

Conversely, negative phylogenetic autocorrelation corresponds to dissimilarities among tips localised in specific parts of the tree, which we call local structures. A local structure would be observed whenever closely related taxa tend to be more different with respect to a given trait than randomly chosen taxa (e.g., Figure 1E). Local structures correspond to relatively recent events that induced divergence of the evolutionary strategies close to the tips of the phylogenetic tree, such as convergence and character displacement (following past or present biotic interactions). This also occurs when the trait under study has been selected towards different optimal values, resulting in opposed evolutionary strategies being observed in sister taxa.

Both structures can be identified using Moran's index (equation 1). The sign of Ecological and life-history strategies of species require not one, but several traits to be adequately described. Accordingly, ecological and life-history strategies are likely to involve combinations of traits with both global and local phylogenetic structures. In the following, we describe a methodology which explicitly investigates multivariate phylogenetic structures that have barely been considered so far.

I

The phylogenetic principal component analysis

Dray et al. (2008) and [START_REF] Jombart | Revealing cryptic spatial patterns in genetic variability by a new multivariate method[END_REF] developed a multivariate approach for identifying spatial structures in multivariate data. Essentially, this approach consists in constraining the principal components of a multivariate method to exhibit spatial autocorrelation, as measured by Moran's index. This methodology proved better at detecting autocorrelated patterns than usual multivariate methods such as principal component analysis [START_REF] Dray | Spatial ordination of vegetation data using a generalization of Wartenberg's multivariate spatial correlation[END_REF][START_REF] Jombart | Revealing cryptic spatial patterns in genetic variability by a new multivariate method[END_REF]. Here, we use the same rationale to define the phylogenetic principal component analysis (pPCA), a method designed to summarise a set of traits into a few synthetic variables exhibiting global or local phylogenetic structures. Note that while we presented pPCA for the analysis of quantitative traits for the sake of simplicity, this approach can be extended to qualitative traits, or even to mixtures of quantitative and qualitative variables [START_REF] Dray | Spatial ordination of vegetation data using a generalization of Wartenberg's multivariate spatial correlation[END_REF].

We denote translates into finding the appropriate loadings u ∈ R p (with u 2 = 1) that minimise and maximise, respectively, the function:

X = [x ij ] (X ∈ R n×p ) a
f : R n×p × R n×n × R p -→ R (X, W, u) -→ var(Xu)I W (Xu) (4)
The solution to this problem is given by the diagonalisation of the matrix 1 2n X T (W + W T )X [START_REF] Dray | Spatial ordination of vegetation data using a generalization of Wartenberg's multivariate spatial correlation[END_REF][START_REF] Jombart | Revealing cryptic spatial patterns in genetic variability by a new multivariate method[END_REF]. It results in a set of loadings {u 1 , . . . , u k , . . . , u r } with u k ∈ R p forming linear combinations of traits (Xu k , the so-called principal components) associated with decreasing eigenvalues λ k , so that:

var(Xu k )I W (Xu k ) = λ k (5)
The largest eigenvalues likely correspond to a large variance and a strong posi- in each principal component (Xu k ) can be computed for a better interpretation of each structure. Moreover, the loadings u k can be used to assess how traits contribute to a given principal component, and thus understand the nature of the corresponding biological structure.

One important choice is that of the phylogenetic weights (W) used in the analysis. Here, we use the measure of phylogenetic proximity underlying the test of [START_REF] Abouheif | A method for testing the assumption of phylogenetic independence in comparative data[END_REF] to define W, because of its good performances at detecting phylogenetic structures [START_REF] Pavoine | Testing for phylogenetic signal in life history variable: Abouheif's test revisited[END_REF]. The phylogenetic proximity a ij among tips i and j is defined as:

a ij = 1 p∈Pij dd p for i = j (6)

A c c e p t e d m a n u s c r i p t

where P ij is the set of internal nodes on the shortest path from tips i to j, and dd p is the number of direct descendants from the internal node p. The phylogenetic proximity a ij defines the entries of the off-diagonal terms of W, all diagonal entries being set to 0. As W is row-standardised, it is defined as:

w ij = a ij n j=1,i =j a ij (7)

Sensitivity study

Extensive simulations were conducted to evaluate the sensitivity of the pPCA to various parameters. Datasets were simulated with different characteristics concerning the type of tree, the tree size, the type, strength and numbers of phylogenetically structured traits, and the total number of traits (including structured and unstructured traits). These parameters are summarised in Table 1. Five types of trees of 16, 32, or 128 tips were simulated to encompass a wide range of tree topologies and sizes: completely symmetric trees (Figure 1A), trees obtained by random clustering of tips (as implemented by rtree function of the ape package, [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF], Figure 1B), the Yule model (Yule, 1924, Figure 1C), the biased model [START_REF] Kirkpatrick | Searching for evolutionary patterns in the shape of a phylogenetic tree[END_REF], Figure 1D), and completely asymmetric trees (Figure 1E). Datasets including random traits and phylogenetically structured traits (i.e., displaying global and/or local structures) were obtained for each tree. Random traits were drawn from a normal distribution (N (0, 1)), while 'structured traits' were obtained by adding normally-distributed random noise to phylogenetic eigenvectors of W [START_REF] Peres-Neto | A unified strategy for estimating and controlling spatial, temporal and phylogenetic autocorrelation in ecological models[END_REF]. Whenever several structures of the same type (global or local)

were simulated in a given dataset, these were derived from the same eigenvector, so that we could evaluate the performance of pPCA when a 'consensus' phylogenetic signal exists in a set of traits (e.g., Figure 1B). This was consistent with the fact that several phylogenetically structured traits are expected to exhibit the same patterns, either because these structures are caused by the same evolutionary process, or because all traits are correlated to another phylogenetically structured trait. 200 datasets were simulated for each of the 810

A c c e p t e d m a n u s c r i p t

different combinations of parameters, resulting in a total of 162 000 datasets.

All simulations were performed in the R software (R Development Core

Team, 2009). Trees were simulated using the R packages ape [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF] and apTreeshape [START_REF] Bortolussi | apTreeshape: statistical analysis of phylogenetic tree shape[END_REF], and scripts developed by TJ for the symmetric model. Structured traits were simulated using the ade4 package [START_REF] Chessel | The ade4 package-I-one-table methods[END_REF][START_REF] Dray | The ade4 package -II: Two-table and K-table methods[END_REF], and data were handled using the phylobase package [START_REF] Bolker | phylobase: base package for phylogenetic structures and comparative data[END_REF].

Each dataset was analysed by a pPCA using Equation 7to define phylogenetic proximities. In each analysis, the structured traits were compared to the Variations in |ρ| according to the different simulation parameters were investigated using a linear model. The relationship between |ρ| and the predictors was linearised using a logit link, i.e. using logit(|ρ|) = log |ρ| 1-|ρ| as the response variable. When interpreting coefficients of the model, predictions μ were re-transposed onto the |ρ| scale, that is, replacing μ by 1 1+e -μ . Qualitative variables were modeled using treatment-coded contrasts (Faraway, 2004, p. 173). The explanatory variables were the type of tree (factor 'tree', the biased model being the intercept), the type of structuring (factor 'strutype', with level 'global' at the intercept), the number of tips ('ntips', intercept=16), the total number of traits ('ntraits', intercept=10), the standard deviation of the random noise added to structured variables ('noise', intercept=0.5), and the number of

A c c e p t e d m a n u s c r i p t

structures (factor 'nstruc', intercept=1).

Empirical data analysis

We illustrate pPCA with the data on reproductive and morphometric traits within a small group of species and subspecies from the lizard family Lacertidae published by [START_REF] Bauwens | Covariation of life-history traits in lacertid lizards: a comparative study[END_REF]. These data are currently available as the dataset lizards in the adephylo package [START_REF] Jombart | adephylo: exploratory analyses for the phylogenetic comparative method[END_REF] for R (R Development Core The analyses were conducted in the R software (R Development Core Team, 2009), using the ade4 package [START_REF] Chessel | The ade4 package-I-one-table methods[END_REF][START_REF] Dray | The ade4 package -II: Two-table and K-table methods[END_REF] for factorial analyses and adephylo [START_REF] Jombart | adephylo: exploratory analyses for the phylogenetic comparative method[END_REF] to perform the pPCA. As in Bauwens and Diaz-Uriarte (1997), data were log-transformed and regressed onto mean adult female length to partial out the body size effect. As a consequence, mean adult female length was removed from the analysis. We then investigated phylogenetic structures of the transformed life-history traits using pPCA.

Results

Sensitivity study

All explanatory variables had a very significant effect on the response variable (Appendix A, Table A , 1997). For a given number of structured traits, the number of random traits slightly lowered the method's ability to retrieve phylogenetic patterns (Δ | ρ| = 0.10 between 10 and 50 traits). Phylogenetic structures incorporating larger amounts of random noise were also more difficult to retrieve (Δ | ρ| = 0.10 between noise of 0.5 and 1). Lastly, the number of structured traits and the type of tree only marginally affected the ability of pPCA to identify phylogenetic patterns.

Empirical data analysis

Both global and local phylogenetic structures were found by pPCA in the lacertid lizards data (Figure 2). The first global principal component of pPCA first opposed a lineage with three species (Lacerta schreiberi, L. agilis, L. vivipara)

having the largest negative scores to the rest of the tree (taxa with positive scores, or scores closer to zero), with the subspecies Podarci h. h. exhibiting the most opposite life histories (Figure 2A). Among the remaining species or subspecies, Lacerta monticola cantabrica and L. m. cyreni were distinctive by their negative scores. The loadings of the analysis (Figure 2B) provided further insights on the corresponding evolutionary strategies, and showed a trade-off between the frequency of reproductive events per year (clutch.F) and the clutch size (clutch.S). L. schreiberi, L. agilis, L. vivipara, and to a lesser extent L. m. cantabrica and L. m. cyreni, reproduce less often but produce a larger number

A c c e p t e d m a n u s c r i p t

of eggs per reproductive event than other populations. A possible explanation for this structure is that environmental conditions only allow for a few reproduction events in these populations, which then deliver lots of eggs with poor individual survival rates. In contrast, the first local principal component of pPCA (Figure 2A) highlighted a strong opposition between related taxa, especially between Takydromus tachydromoides and Acanthodactylus erythrurus, but also to a lesser extent between L. schreiberi and L. agilis.

Discussion

Phylogenetic autocorrelation has so far been considered as a mere nuisance to the correlative analysis of comparative biological data, when exploring tradeoffs as well as in allometric studies. In this paper, we advocate that phylogenetic autocorrelation is a source of relevant biological information for the exploratory analysis of such data. To accomplish this task, we introduced the phylogenetic principal component analysis (pPCA), a method that we adapted from existing multivariate spatial statistics [START_REF] Dray | Spatial ordination of vegetation data using a generalization of Wartenberg's multivariate spatial correlation[END_REF][START_REF] Jombart | Revealing cryptic spatial patterns in genetic variability by a new multivariate method[END_REF] to analyse phylogenetic structures in multivariate sets of traits. Based on the results obtained from simulated and empirical data, we discuss the ability of the pPCA to retrieve phylogenetic signals in a multivariate set of traits, and the impact that this approach could have in evolutionary ecology. Therefore, local structures may be interpreted even though the corresponding eigenvalue seems negligible compared to global structures, provided it is biologically significant. Other results of our sensitivity study suggest that pPCA performs better in larger trees, although performances on small phylogenies were satisfying. Interestingly, pPCA seemed rather insensitive to the shape of the, indicating that the method can be used with virtually any kind of phylogeny.

Although pPCA will be best appreciated using empirical datasets, further simulation studies may be considered. In this study, we used eigenvectors of a phylogenetic proximity matrix [START_REF] Peres-Neto | A unified strategy for estimating and controlling spatial, temporal and phylogenetic autocorrelation in ecological models[END_REF] to simulate phylogenetic structures. This method allows simulation of complex phylogenetic patterns in negligible computational time, which permitted examination of the influence a large number of parameters on pPCA results. The drawback of this approach is that eigenvectors of phylogenetic proximity matrices are not directly related to an model of evolution such as the Brownian motion or the OU models. While current procedures implementing trait simulation under these models are more computer-intensive, it could be possible to study how pPCA behaves under these models, given variation in a few parameters.

The main parameter that should be investigated in further detail is the

A c c e p t e d m a n u s c r i p t

phylogenetic proximity used in pPCA. A previous study demonstrated that some phylogenetic proximities were better than others at detecting phylogenetic structures using Moran's index [START_REF] Pavoine | Testing for phylogenetic signal in life history variable: Abouheif's test revisited[END_REF]. However, it is likely that the most appropriate measurement of phylogenetic proximity depends on the dataset under scrutiny. To assess whether a given phylogenetic proximity is adapted to a particluar dataset, we advocate to perform Moran's I test using this proximity matrix. Whenever significant structures are detected, one can input this phylogenetic proximity in pPCA to uncover the nature of the underlying phylogenetic structures.

Potential impacts of the approach in evolutionary and ecological studies

This novel approach should complement nicely the usual PCM toolbox, Life histories, for example, are likely to be phylogenetically structured [START_REF] Gaillard | An analysis of demographic tactics in birds and mammals[END_REF][START_REF] Pontier | Maternal investment per offspring and demographc tactics in placental mammals[END_REF][START_REF] Rochet | Comparative analysis of phylogenetic and fishing effects in life history patterns of teleost fishes[END_REF]. In our case study, the pPCA identified phylogenetic patterns in the main life-history tactics adopted by a set of taxa. Our results suggest that the trade-off between clutch frequency and size may have resulted in the ancient divergence of evolutionary strategies. In contrast, the trade-off between hatchling mass on the one hand as speciation caused by diversifying selection or niche separation, and are thus also of fundamental interest. Dating these recent events would be even more interesting as historical information about the considered taxa is more likely to be available for recent speciation events. For instance, we could investigate whether a recent speciation highlighted as local structure would have been preceded by significant modifications of the environment.

A further strength of pPCA lies in its ability to analyse very large sets of traits (i.e., hundreds or thousands of traits) simultaneously. Usual PCM typically rely on pairwise comparisons among traits, which becomes cumbersome when lots of variables are under scrutiny, and often requires discarding traits from the analysis. This issue will be increasingly concerning in the near future as new and large databases of life-history traits will become available. pPCA can be used to explore such data, to unveil evolutionary trade-offs among a large number of traits, without having to make a prior selection of analysed traits.

Previous methods have already attempted to analyse phylogenetic signals in a series of traits. These methods determine the proportion of variation in a set of traits correlated with the phylogenetic relatedness among species (Giannini,

A c c e p t e d m a n u s c r i p t

2003; [START_REF] Desdevises | Quantifying phylogenetically structured environmental variation[END_REF]. They involve factorial analyses, but with the aim of partitioning variation in traits instead of depicting phylogenetic structures.

Nevertheless, one step of the [START_REF] Giannini | Canonical phylogenetic ordination[END_REF] variation partitioning approach consists of selecting the nodes of the phylogeny that better explain the variation in the trait values. Accordingly, the selected nodes could be used to determine at which depth in the phylogeny the taxa differ in their trait values. However, this selection results from a long series of tests to determine if the differences between the lineages that descend from a given node are significantly responsible for trait variation. The number of tests depends on the size of the phylogenetic tree, with an increasing risk of erroneously significant tests. The pPCA approach thus brings a new optimised way of disentangling the phylogenetic patterns in a set of traits by identifying the lineages and also the combination of traits responsible for global versus local trait variation.

A potential application of the pPCA concerns phylogenetic community ecology [START_REF] Webb | Phylogenies and community ecology[END_REF]. Phylogenetic clustering in a community (lower phylogenetic diversity than expected by chance in a regional pool of species) merely reflects the simultaneous action of environmental filtering and phylogenetic conservatism. In contrast, distinct, even opposed processes can lead to phylogenetic overdispersion (higher phylogenetic diversity within a community than expected by chance in a regional pool of species). For example, phylogenetically overdispersed communities can arise from (i) limiting similarity and conservative traits or (ii) environmental filters with convergent traits [START_REF] Kraft | Trait evolution, community assembly, and the phylogenetic structure of ecological communities[END_REF]. Consequently, knowledge of the evolution of traits is necessary to interpret observed structures in phylogenetic diversity. A difficulty is that the traits involved in environmental filters and those involved in limiting similarity might follow different evolutionary pathways [START_REF] Emerson | Phylogenetic analysis of community assembly and structure over space and time[END_REF][START_REF] Ackerly | Niche evolution and adaptive radiation: testing the order of trait divergence[END_REF].

Although pPCA does not provide a formal test of phylogenetic conservatism or over-dispersion, it can be used to describe the phylogenetic signal induced by these processes. Therefore, our methodology could be applied to describe the level of phylogenetic signal in sets of traits from labile traits with local close-

A c c e p t e d m a n u s c r i p t

to-tips variations to conserved traits with global close-to-root variations. This approach could provide even more details by highlighting lineage-dependent signals. For example, a single trait might exhibit a general pattern of phylogenetic signal (global structure) and also strong localised trait variations in a single lineage (local structure). Moreover, pPCA can be used together with other multivariate methods to relate the combination of traits identified with a given phylogenetic structure (either global or local) to explanatory factors, such as environmental variables. For instance, co-inertia analysis [START_REF] Dolédec | Co-inertia analysis: an alternative method for studying species-environment relationships[END_REF]Dray et al., 2003a,b) could be used to link phylogenetic structures identified by pPCA to descriptors of the ecological niche, so as to assess potential patterns of adaptation. pPCA could therefore complement both existing ecological methods (e.g., co-inertia) and evolutionary approaches (e.g., phylogenetic overdispersion/clustering), providing a link between trait evolution, patterns in phylogenetic diversity, and biotic or abiotic interactions, and giving insights into the historical and ecological processes that underpin community assembles.

To conclude, we illustrate the intersection between issues in spatial and phylogenetic methods. Spatial and phylogenetic patterns are generated by very different processes, but the mathematical tools that can be used to measure and model these patterns may be similar. This is because both rely on the concept of autocorrelation, which can be defined as the non-independence among observations with respect to a set of underlying proximities. Several spatial methods developed in ecology have already been successfully adapted to PCM (Cheverud et al., 1985;[START_REF] Gittleman | Adaptation: statistics and a null model for estimating phylogenetic effects[END_REF][START_REF] Diniz-Filho | An eigenvector method for estimating phylogenetic inertia[END_REF][START_REF] Desdevises | Quantifying phylogenetically structured environmental variation[END_REF][START_REF] Giannini | Canonical phylogenetic ordination[END_REF]. Originally, spatial autocorrelation was perceived by ecologists as a nuisance that precluded the use of standard statistical tools in correlative studies. However, the study of spatially autocorrelated patterns turned out to be a fecund paradigm, as ecologists realised that these structures were mere indicators of considerable underlying ecological processes.

The same may be true of phylogenetically autocorrelated patterns. Rewording [START_REF] Legendre | Spatial autocorrelation: trouble or new paradigm?[END_REF], we can now also ask the question: is phylogenetic autocorrelation trouble, or a new paradigm?

A c c e p t e d m a n u s c r i p t

Acknowledgements

We are grateful to the CCIN2P3 for providing access to their computers, and particularly to Simon Penel for his help. We thank Anne-Béatrice Dufour for useful discussions on a former version of the manuscript. We address many thanks to F. Stephen Dobson and Nigel G. Yoccoz for their insightful comments and thorough review of our manuscript. We finally wish to thank the team developing the phylobase package for easing the handling of phylogenetic comparative data. (2) number of phylogenetically structured traits (global/local).

(3) standard deviation of normal variates added to phylogenetically structured traits. (4) total number of traits in the dataset, including phylogenetically structured traits. standard deviation of the random noise added to the structured traits.

  Phylogenetic autocorrelation relates to the non-independence of trait values observed in taxa given their phylogenetic proximity. There are two ways in which this non-independence can arise, depending on whether closely related taxa tend to have more similar, or more dissimilar trait values than expected at random, resulting in positive and negative autocorrelation, respectively. Positive phylogenetic autocorrelation most often results in global patterns of similarity in related taxa; we thus refer to these patterns as global structures. Global patterns reflect the general idea of phylogenetic signal: trait values observed in a set

  depends on how values of a trait (x i ) relate to the values observed on closely related taxa (Wx i ). Moran's I will be greater than (respectively less than) I 0 (value of I in the absence of autocorrelation) when closely related taxa tend to have similar (respectively dissimilar) values for the studied trait. Obviously, the definition of phylogenetic proximities in W will condition the measurement of global and local structures. As shown by[START_REF] Pavoine | Testing for phylogenetic signal in life history variable: Abouheif's test revisited[END_REF], not all phylogenetic proximities are equal in detecting phylogenetic structures. Especially, the phylogenetic proximities underlying Abouheif's test (matrix A = [a ij ] in[START_REF] Pavoine | Testing for phylogenetic signal in life history variable: Abouheif's test revisited[END_REF]) proved superior to several common phylogenetic proximities for testing phylogenetic inertia in traits simulated under Brownian and Ornstein-Uhlenbeck (OU) processes. More generally, the matrix W can be derived from any model of evolution which seems appropriate to the data, taking branch lengths into account whenever these are accurately estimated, and rely-the topology in other cases.

  matrix containing p quantitative traits measured on n taxa, and W a matrix of phylogenetic weights used in the computation of Moran's I (equation 1). As in classical PCA, missing data can be set to the mean of the corresponding trait, which does not add artefactual structures to the analyzed traits. Without loss of generality, we assume that traits are centered (i.e., i x ij = 0 with j = 1, . . . , p). The purpose of pPCA is to find linear combinations of traits (columns of X) containing a large variance and displaying global or local phylogenetic structures. Mathematically, this problem

  tive I, indicating global structures (close-to-root variation in trait states). Conversely, the lowest (i.e., most negative) eigenvalues correspond to a high variance and a large negative I, indicating local structures (close-to-tips variation in trait states). As in other reduced space ordinations, the eigenvalues indicate the amount of structure expressed by each synthetic variable. A sharp decrease in the screeplot is likely to indicate a shift between strong and weak structures. The amount of variance (var(Xu k )) and phylogenetic autocorrelation (I(Xu k ))

  first relevant (global and/or local) principal component of pPCA, to assess how the method performed. The strength of the link between the original simulated structures and patterns identified by pPCA was measured using the absolute value of Spearman's rank correlation, |ρ|. Whenever the dataset included several distinct structured traits, |ρ| values were averaged by type of structure (i.e., global or local). Hence, we obtained one or two |ρ| per simulated dataset, used as indicator of the performance of pPCA (|ρ| close to one = high performance, |ρ| close to zero = low performance).

  Performance of the approach -methodological discussionPreliminary results stemming from our sensitivity study were very promising, and provided guidelines for applications of pPCA. Overall, pPCA performed well to retrieve simulated phylogenetic structures, even in some cases where only 1 out of 50 traits was phylogenetically structured. pPCA seemed to retrieve global phylogenetic structures more easily than local structures. This may be due to the asymmetry of Moran's I distribution, which often has a smaller range of variation in negative values (local structures) than in positive values (global structures)[START_REF] De Jong | On extreme values of Moran's I and Geary's c[END_REF]. As pPCA seeks principal components with extreme values of I, global structures (associated with large positive I) would be more easily detected than local structures (associated with large negative I).

  bringing a new perspective to the analysis of comparative biological data. Contrary to usual PCM, our approach does not attempt to improve estimates of correlations among traits by 'correcting' for phylogenetic dependence among species. Instead, it seeks biologically meaningful combinations of traits that are globally or locally phylogenetically structured, thus allowing us to uncover fundamental evolutionary patterns. As noted by Bauwens and Diaz-Uriarte (1997), theories of life-history evolution are explicitely micro-evolutionary [...] whereas patterns of life-history covariation are most evident when comparisons are made among higher taxonomic levels. pPCA covers both of these aspects, by providing insights about broad macro-evolutionary patterns (global structures) and more recent, even micro-evolutionary patterns (local structures).

  and frequency on the other hand appears to be more labile, involving more recent character changes in most of the lineages and especially between T. tachydromoides and A. erythrurus. The pPCA thus allows description of global (close-to-root, phylogenetic signal) versus local (close-to-tips) phylogenetic structures in a multivariate set of traits, and highlights which lineages and which taxa are involved in these structures. Overall, this illustration using an empirical dataset showed that pPCA can bring new insights about evolutionary strategies of a set of taxa. Moreover, whenever a molecular clock is available for the considered phylogeny, it would be possible to estimate the age of the involved lineages and taxa, by dating their most recent common ancestors. This would allow assessing how and when different evolutionary strategies might have appeared and evolved along the history of the considered taxa. Local structures uncovered by pPCA point out more recent evolutionary events, such

Figure 1 :

 1 Figure 1: pPCA of simulated data. Example of simulated traits (light yellow) for different tree structures (A-E), and structures identified by pPCA. Global components and corresponding eigenvalues are indicated in red, while local components and their eigenvalues are displayed in blue. Positive and negative values of traits and PCs are indicated by black and white circles, respectively. Symbol size is proportional to absolute values. Simulated traits are labelled as Gi: i th global structure, Li: i th local structure, and Ri: i th random (i.e., non-phylogenetically structured) trait. Principal components (PC) of pPCA are labelled as GPC1: first global PC (i.e., associated with the largest positive eigenvalue). LPC1: first local PC (i.e., associated with the largest negative eigenvalue). (A) Symmetric tree; random noise added structures ('noise') equaled 0.5. (B) Random clustering of tips; noise=1. (C) Yule model; noise=0.5. (D) biased model; noise=0.75. (E) Assymetric tree; noise=1.

Figure 2 :

 2 Figure 2: pPCA of lizards data. (A) First global (red) and local (blue) principal components of the pPCA of lacertid lizards data, after removal of size effect. Inset barplot displays the corresponding eigenvalues. Positive and negative scores are indicated by black and white circles, respectively. Symbol size is proportional to absolute values. Taxa are labelled as: Podarcis h. atrata ('Pa'), P. h. hispanica ('Ph'), Lacerta lepida ('Ll'), L. monticola cantabrica ('Lmca'), L. m. cyreni ('Lmcy'), Podarcis h. hispanica Asturias ('Phha'), P. h. h. Salamanca ('Pha'), P. bocagei ('Pb'), P. muralis ('Pm'), Acanthodactylus erythrurus ('Ae'), Takydromus tachydromoides ('Tt'), T. septentrionalis ('Ts'), Lacerta vivipara ('Lviv'), L. agilis ('La'), L. schreiberi ('Ls'), and L. viridis (''Lvir'). (B) Loadings of the traits for the first global (red) and local (blue) principal components. Inset barplot displays the corresponding eigenvalues. d=0.5 indicates the mesh of the grid. Analysed traits are hatchling length (hatch.L) and mass (hatch.m), clutch frequency (clutch.F) and size (clutch.S), mean and

  Figure 2:

( 1 )

 1 trees obtained by random clustering of tips. (2) comb-like model (completely asymmetric trees). (3) completely symmetric trees. (4) Yule model. (5) local phylogenetic structure. (6) number of tips. (7) total number of traits. (8) number of structured traits (1 or 3). (9)

Table legends Table 1 :

 legends1 parameters of the simulated data. 200 datasets were simulated for all combinations of these parameters. (1) expressed in number of tips.

  Table A.2: Coefficients of the model. Factors are preceded by 'fac', followed by the levels.

		Estimate Std. Error t value Pr(>|t|)
	(Intercept)	0.8474	0.0053	158.99 < 2e -16
	fac.treeclust 1	0.0431	0.0053	8.16 < 2e -16
	fac.treecomb 2	0.2980	0.0053	56.48 < 2e -16
	fac.treesym 3	0.4450	0.0053	84.34 < 2e -16
	fac.treeyule 4	0.0410	0.0053	7.78 < 2e -16
	fac.strutypelocal 5	-1.4759	0.0033 -442.31 < 2e -16
	ntips 6	0.0046	0.0000	137.77 < 2e -16
	ntraits 7	-0.0223	0.0001 -227.45 < 2e -16
	fac.nstruc3 8	0.0845	0.0033	25.32 < 2e -16
	noise 9	-0.9520	0.0082 -116.47 < 2e -16

(2) type of structure (global or local). ( 3) number of tips. (4) total number of traits. ( 5) number of structured traits (1 or 3). ( 6) standard deviation of the random noise added to the structured traits.