How do functional traits syndromes covary with growth and reproductive performance in a water-stressed population of Fagus sylvatica?

A central issue in plant evolutionary ecology is to understand how several coordinated suites of traits (i.e. traits syndrome) may be jointly selected within a single species. This study aims to describe patterns of variation and co-variation of functional traits in a water-stressed tree population and test their relationships with performance traits. Within a Mediterranean population of Fagus sylvatica experiencing recurrent summer droughts, we investigated the phenotypic variation of leaf unfolding phenology, leaf area (LA), leaf mass per area (LMA), leaf water content (LWC), water use efficiency (WUE) estimated by carbon isotopic discrimination (d13C), twig Huber-value (HV: the stem cross-section divided by the leaf area distal to the stem), wood density (WDens), and leaf nitrogen content (N mass ). First, a principal component analysis revealed that two main axes structured the phenotypic variability: the first axis opposed leaf unfolding earliness and LWC to LMA and WUE; the second axis opposed LA to HV. These two axes can be interpreted as the opposition of two strategies (water economy versus water uptake) at two distinct scales (leaf for the first axis and branches for the second axis). Second, we found that LMA, LA, leaf unfolding and LWC responded differently to competition intensity, while WUE, WDens and HV did not correlate with competition. Third, we found that all studied functional traits were related to growth and/or reproductive performance traits and that these relationships were frequently non-linear, showing strong interactions between traits. By


Introduction
A major issue of plant evolutionary ecology is to understand how within-species phenotypic variation of functional traits allows the adaptation of populations to their environment.Populations can adjust their traits to environmental variation both through a plastic response within the course of individual life and/or through a genetic response to selection across generations.Accordingly, these two responses contribute to shape the standing phenotypic variation at a given moment in time.The direct study of phenotypic variation of functional traits in nature and their relation to fitness can bring valuable information on adaptive processes (Dudley 1996, Heshel andRiginos 2005).Several morphological, physiological or phenological traits simultaneously affect fitness via their effects on growth, survival and fecundity (Violle et al. 2007).Rather than single traits, natural selection is thus likely to target coordinated suites of traits (often called syndromes), and to favour one or several ecological strategies, i.e. sets of individual trait values that increase fitness within a given environment (Arntz and Delph 2001, Donovan et al. 2011, Laughlin and Messier 2015).
The study of traits syndromes, of their contribution to adaptation and of their consistency within and among species is a growing subject of interest for ecologists (Albert et al. 2010, Wright andSutton-Grier 2012).When studying the multivariate pattern of variation of functional traits within a single environment, the syndromes associated to each axis of variation define alternative strategies.Interspecific studies are often used to map multivariate correlation patterns of functional traits (e.g.Wright et al. 2004), considering that different species sharing the same environment represent different adaptive solutions to ecological conundrums (Wright et al. 2004, Shipley et al. 2006, Díaz et al. 2016).On the other hand, there is increasing evidence that functional trait variation is very important within species as well as within populations (Alberto et al. 2013) and thus investigating how traits cluster and affect plant performance at these biological scales can provide important insights on how strategies are selected (Donovan et al. 2011).
Regarding leaf morphological and physiological traits, two opposed strategies, resource acquisition versus resource conservation, were found to explain a large part of the observed syndromes in interspecific studies (Wright et al. 2004).In the case of adaptation to seasonal water deficit, two main ecological strategies are known: (1) a water economy strategy, where plants maintain low growth rates and low rates of gas exchange during droughts, and (2) a water uptake strategy, where plants have a more rapid instant growth through higher rates of gas exchange when water is available, typically spring in Mediterranean climate, allowing them to complete important biological functions before drought onset (Lo Gullo and Salleo 1988, Ludlow 1989, Arntz and Delph 2001).The phenology of leaves development is also likely to interact with leaf morphological and physiological traits, and to play a role in the water fluxes (i.e.longer leafy period increases transpiration).Within-species however, and in particular for long-lived species within a drought-prone environment, it is unclear whether these strategies are meaningful, and whether selection could favour intermediate phenotypes and/or different trait combinations.
To investigate the adaptive significance of traits syndromes, one solution is to relate individual variation of functional traits to fitness within a population (Donovan et al. 2011).The classical approach is the phenotypic selection analysis that allows to model the so-called fitness surface (Dudley 1996, Heshel andRiginos 2005) by characterising the relationship between traits and fitness using univariate or multivariate linear

Version postprint
Comment citer ce document : Bontemps, A., Davi, H., Lefèvre, F., Rozenberg, P., Oddou-Muratorio, S. (2017).How do functional traits syndromes covary with growth and reproductive performance in a water-stressed population of Fagus sylvatica? .Oikos., DOI : 10.1111/oik.04156 'This article is protected by copyright.All rights reserved.' regression models (e.g.Lande and Arnold 1983).In both the univariate and multivariate models, the selection on a single trait X is classically summarised by two parameters: the directional and the quadratic selection gradients (β X and γ X ), that respectively characterize the slope and shape of the fitness surface.In the multivariate model, a third parameter, the correlational selection gradient (β XY ), can then be added to inform on the interactive effect of traits X and Y on fitness.In the linear regression theory, the parameters of the multivariate model are partial regression coefficients measuring the strength of the relationships between traits and fitness while holding constant the other traits (or some other confounding variable), and they can be used to tease apart the direct effect of each traits from their indirect effects through correlated traits.Therefore, the multivariate model is a model of choice to evaluate selection on phenotypically correlated traits (Lande & Arnold 1983).
Recent studies have investigated the variation of tree responses to drought at individualor population-level across vegetation types and environmental conditions (e.g.Choat et al. 2012;Klein et al. 2014, Martínez-Vilalta et al. 2009), but rarely tested their adaptive significance (but see Ramirez-Valiente et al. 2010, Laforest-Lapointe et al. 2014).A major limitation for phenotypic selection studies in trees is that fitness estimation is constrained by their long life cycle.However, some performance traits are recognized as good proxies of survival and reproductive components of fitness in trees.Growth is often used as a predictor for survival (Cailleret et al. 2016).Even if growth correlates with reproduction (Niklas and Enquist 2003), traits such as seed output and seed mass are usually considered as better proxies of reproduction (Clark et al. 2007).
Here, we applied an in-situ approach to infer how functional traits associations contribute to fitness within a Fagus sylvatica L. (European beech) population growing on shallow soils in the southern margin of the distribution range and regularly facing summer droughts in a Mediterranean climate.We studied the variation and co-variation patterns of a set of leaf functional traits and of one wood trait: leaf mass per area (LMA), leaf area (LA), leaf water content (LWC), intrinsic water use efficiency (WUE) estimated by carbon isotopic discrimination (d13C), twig Huber-value (HV), leaf nitrogen content (N mass ) and phenology of leaf unfolding, plus wood density (WDens).
We also characterised several fitness-related performance traits: radial growth rate and fruit output (number and weight of fruits).
In this study of a Fagus sylvatica population in a drought prone environment, our major hypothesis was that antagonistic selective pressures would stem from the maximisation of carbon acquisition versus the capacity to cope with strong summer droughts.For instance, early tree phenology is expected to increase the photosynthetic period during springtime when there is no water stress.But higher photosynthesis during that period could also increase the yearly transpiration and the drought intensity during summer (Davi et al. 2006).Moreover, higher transpiration would allow higher photosynthesis during the period without water stress, at the cost, however, of lower water use efficiency.Besides, other environmental pressures can shape functional traits and may reinforce or counteract adaptation to water stress.For example, late frosts would more likely damage early leafing trees (Menzel et al. 2015) while greater winter embolism would delay leafing (Beikircher et al. 2016).
To investigate this hypothesis of antagonistic selection pressures, we first characterised the multivariate phenotypic correlation structure.Secondly, we tested whether competition intensity (a major environmental factor expected to drive functional traits variation) and spatial autocorrelation contributed to the observed phenotypic correlations.Thirdly, using phenotypic selection analysis, we tested the relationships between functional traits and fitness (measured by performance traits) and investigated directional and antagonistic selection.In particular, we searched for traits synergic effects i.e. whether traits allowing water economy on one side or water uptake on the other side cluster preferentially together or interact to promote performance.Finally, combining these different approaches, we addressed the following questions: 1) What is the phenotypic variance of functional traits in this marginal population, how does it vary spatially and does it correlate to local variation of competition?2) Do multivariate patterns of phenotypic variation reveal the co-existence of several ecologically consistent phenotypic combinations?3) To what extent do the functional traits and their interactions relate to individual fitness?

Study site and species
The study site is located on Mont Ventoux (44°10'28''N; 5°16'16"E, France), where F. sylvatica ranges from ~ 900 m to 1700 m in altitude.The study plot (area: 0.83 ha) is at the lower limit of the local range (987 m to 1048 m).The mother bedrock is calcareous, and the soil is characterised by a highly variable depth and a high proportion of stones (Cailleret et al. 2014) resulting on an average soil water content of 262 mm (with 50% of stone content and 35 cm of soil depth).
Fagus sylvatica is a shade-tolerant species sensitive to water deficit during summer (Lebourgeois et al. 2005, Lendzion andLeuschner 2008).Moreover, F. sylvatica is sensitive to late spring frosts that limit early leaf unfolding (Dittmar and Elling 2006).At the study site, the climate (mountainous Mediterranean) is characterised by summer droughts and mild winters, with rainfalls mostly occurring in autumn (mean annual temperatures: 9.32°C, mean annual precipitations: 1020 mm, mean summer precipitations: 186.80 mm, Supplementary material Appendix 1, Tab.A1).The main tree species expected to compete with F. sylvatica at the study site are Abies alba, Pinus sylvestris, Acer opalus, Ilex aquifolium, Quercus pubescens, Sorbus aria, and Pinus nigra (Fig. 1).Within the site, we exhaustively mapped 175 dominant adult F. sylvatica trees with diameter > 12 cm and with crowns in the upper canopy receiving full sun (Fig. 1).Because vegetative reproduction occurs through suckering or resprouting in F. sylvatica, we checked for the presence of clones using genetic analysis (Bontemps et al. 2016).We found 26 pairs of individuals with the same genotype and only kept the biggest individual within each genet, leading to a total of 149 focal non-clonal dominant trees (Fig. 1).We recorded their spatial coordinates using a Trimble GPS (precision < 1 m) and measured diameter at 1.3 m above ground (dbh) and height.The mean age estimated by counting the number of radial rings at 1.3 m was 92.7 years (s.d.=13.7 years).Finally, we exhaustively mapped all competing individuals in a radius of 20m around each of the 149 focal individuals, and recorded their species and dbh.

Functional traits choice
We measured eight functional traits (Table 1).Their potential role in the response to water deficit has already been highlighted, either through intra-specific studies of plastic variation in response to water deficit (1,3,4,5,6,7,8 below), or through inter-specific studies in water-stressed versus non-stressed habitat (2,3,8 below): Accepted Ar ticle (1) Leaf unfolding determines the vegetation length, therefore the carbon uptake, and the tree transpiration at yearly time scale (Davi et al. 2006).A drought event can delay budburst in the subsequent spring (Kuster et al., 2014) but an earlier budburst could increase transpiration over the current year and consequently amplify summer drought effects (Bréda et al., 2006).(2) A general tendency for species inhabiting arid and semi-arid regions to have high-LMA, leathery leaves has been reported (Wright et al. 2004;Poorter et al. 2009).
(3) Low LA is commonly correlated to low temperature at high elevation for instance (Bresson et al. 2011), but also with high water stress (Ackerly et al. 2002).(4) d13C is considered as a proxy of WUE (Farquhar & Richards 1984) that mechanistically increases when stomata close during a drought.(5) LWC exhibits seasonal variation, with lowest values obtained during summer (Lo Gullo and Salleo 1988).(6) HV increases with drought especially for pines (Delucia et al. 2000;Martínez-Vilalta et al. 2009).(7) Leaf nitrogen content, measured as the quantity of nitrogen per unit of dry leaf mass (N mass ) is related to photosynthetic capacities because RUBISCO represents almost 50% of nitrogen in leaves (Shipley et al. 2006).(8) Finally, wood density correlates with resistance to embolism at interspecific level (Hacke et al. 2001;Chave et al. 2009) and at inter-annual scale, wood density also increases during drought events (Bouriaud et al. 2005).
We summarized the expected sign of phenotypic variation according to water stress gradient for each trait in Table 2.

Sampling methods for functional traits
Leaf unfolding phenology: In 2009, we recorded in-situ leaf unfolding dynamics.Five stages were distinguished, noted from 1 to 5 (Davi et al. 2011).The phenological stages of the upper and lower part of each adult tree were noted on 15 different dates from the 23 rd of March to the 4 th of May, and then averaged at tree level.A single Phenology Score Sum (PSS) was computed for each tree as the sum of the phenological stages observed over all of the dates.The higher the PSS at a given date of measurement, the earlier and quicker was leaf unfolding.

Morphological and physiological traits:
In the second week of July 2009, F. sylvatica trees were climbed to collect light-exposed branches at the top of the crown.The branch part carrying the last growth unit and ten terminal leaves was cut.For each individual, three branches were sampled at different locations within the top of the crown and treated as replicates.All sampled branches were directly placed in plastic bags within a cooler in the field and were taken to a cold chamber (+4°C) at the end of the sampling day.
Global fresh mass and area of the ten terminal leaves were measured for each of the three replicates.Leaf area (LA) was measured using a leaf area meter (Delta-T Area meter, Delta-T Devices Ltd., Cambridge, UK).Then, leaves were oven-dried at 65°C during 48 h and global dry mass (LM) was measured.Leaf Water Content (LWC) was calculated as the ratio of fresh mass minus dry mass over fresh mass and leaf mass per area (LMA) as the ratio of dry mass over leaf projected area.We then assessed the

Version postprint
Comment citer ce document : Bontemps, A., Davi, H., Lefèvre, F., Rozenberg, P., Oddou-Muratorio, S. ( 2017).How do functional traits syndromes covary with growth and reproductive performance in a water-stressed population of Fagus sylvatica? .Oikos., DOI : 10.1111/oik.04156 'This article is protected by copyright.All rights reserved.' Huber Value (HV) as the stem cross-section divided by the leaf area distal to the stem (based on the last growth unit).Among the three replicates per tree, the one with the highest LMA (expected to be the most light-exposed, Davi et al. 2009) was selected for carbon isotope and nitrogen analysis (3 leaves mixed together).Carbon isotope composition ( 13 C and 12 C) of each sample was measured using a continuous flow isotope mass spectrometer coupled to an elemental analyser (Roussel et al. 2009).Sun leaves used CO 2 from the top of the canopy, where the atmosphere is well mixed and therefore we assume that all trees received equivalent CO 2 concentrations and that d13C values are therefore comparable among individuals.
Trees were cored at breast height at the end of the winter of 2008-2009.The samples were dried to 12% water content and X-rayed.The X-ray images were scanned at 4000 dpi and wood density (WDens) was computed as the average density over the whole core.

Performance traits
We measured four performance traits considered as components of fitness (Table 1).

Growth:
The mean ring width over the last ten years (period 1998-2008), noted as RW10M, was used to characterise growth.For 98 individuals, we measured ring width by using wood cores X-ray images and WinDENDRO software (Guay et al. 1992).For the other 51 individuals, the wood cores X-ray images were blurred (due to strong grain angle); however we could distinguish the rings using scanned images and CDendro v5.3 software (Larsson 2003).We checked for the presence of missing rings and pointing errors by visually cross-dating the individual ring width profiles with a reference chronology available from Cailleret et al. (2014).We kept 137 individual profiles out of 149 for ring width analysis.The RW10 variable was calculated using the last ten rings of each tree.These rings are weakly affected by low-frequency internal and external factors like cambial age and long-term climate variation.Hence detrending techniques used in dendrochronology for long time-series were not relevant in our case.
Reproductive traits: Seed production was characterized by relative fruit output at tree level (FO tree ), fruit output at branch level (FO branch ) and mean fruit weight (FW mean ).We estimated visually the total fruit output at the end of summer 2009 (masting year) through the percentage of the crown covered by fruit shells according to the following classes: 1) no fruits observed, 2) less than 25%, 3) between 25% and 50%, 4) between 50% and 75%, 5) over 75%.The fruit output per branch and fruit weight were assessed on the same branches collected for leaf trait analysis, with a mean number of fruits per branch of 2.1 fruits and a standard deviation equal to 1.5 fruits.

Statistical analyses
Correlation between functional traits: we used Principal Component Analysis (PCA) to characterise the association among the eight studied functional traits (Table 1) at individual scale and identify the trait syndromes.Traits were scaled to mean zero and unit variance prior to the analyses.As eight traits were included in the PCA, we considered only PCA axes capturing at least 1/8 of the total variation (here at least 12%).Then two supplementary variables were projected on the principal components (PCs): tree height (Height) and the competition index estimated at the distance interval of 20 m (Comp20), in order to control how both variables covary with the other functional traits.The analysis was performed with the function PCA of the R software package FactoMineR (Lê et al. 2008).Pearson coefficients of correlations were also computed and tested.Effect of competition on functional traits: We used the Martin-Ek index (Martin and Ek 1984) to quantify the intensity of intra-specific and inter-specific competition on a focal individual i.This index accounts simultaneously for the diameter (dbh) and the distance of each competitor j to the competed individual i:

Accepted Ar ticle
Eq. 1 where d bhi and d bhj are the diameter at breast height (in cm) of the competed individual i and of competitor j (any adult tree of any species with dbh j >dbh i ), n dmax the total number of competitors in a given radius d max (in m) around each individual i , and d ij the distance between individuals i and j.Overall, we considered competitors in a radius from 5 to 20 m (Comp 5 to Comp 20 ).Note that this index assumes that intra-and interspecific competition are equivalent.Spatial autocorrelation of traits: we investigated the relationship between phenotypic si ilarity and physical distance between pairs of individuals using the Moran's I statistics and the software AUTOCORQ (Hardy 2009).These spatial auto-correlation analyses were done for each functional and performance trait, and also for the four first P 's identified by the P A. The level of spatial autocorrelation was assessed for each trait by ( ) the regression coefficient of Moran's I statistics on pairwise physical distance (Slope); and (2) the distance up to which Moran's I statistics was significantly positive (D Tresh ).A significant slope indicates that neighbour individuals are significantly more similar phenotypically than distant individuals, up to the DTresh value.Confidence intervals around Slope and D Tresh values were obtained by 5000 random permutations of trait values over individual positions.

Phenotypic selection analyses:
To characterize the selection acting on functional traits, we performed multiple regression analysis of each performance trait against functional traits, and against the four PCs from the PCA described above.For each performance trait, we estimated the relative fitness associated to each performance trait by dividing each individual value by the population mean (Lande & Arnold 1983).Moreover, two covariates were added in these models.First, Comp20 was included to remove the effect of competition on the inferred selection gradients.Second, dbh was included to remove geometric effects on ring width.
For analyses on functional traits, we considered the model: P ∑ ∑ ∑ ∑ i o p 2 dbh , Eq. 2 Where P is the performance trait, α is the origin of the regression, β i is the directional selection gradient on trait X i , γ i is the quadratic selection gradient on trait Xi, β ij are the correlational selection gradients involving interaction effect between traits X i and X j , 1 and 2 the effect of Comp20 and dbh, and the residual.Interactions between competition and traits were also included in the models, in order to account for a potential change in relationships between traits and performance depending on the level of competition.With this complete multivariate model, we aimed to characterise the potential complexity of the selection surface.We also fitted univariate and simpler multivariate regression models on functional traits (results detailed in Supplementary material, Appendix 3).For analyses on PCs, we considered the same model, with trait X i replaced by PC i .These models were fitted using the lm function implemented in R-base (R Core Team, 2016), and following a three-step procedure: 1) We selected the best set of explanatory variables using a stepwise regression procedure, based on the AI (the "step" function of the R software).Prior to this analysis, LA, dbh, RW10M were log-transformed to achieve normality in order to perform efficient parameter selection.Moreover, the functional traits and the competition index were scaled to mean zero and unit variance.Such transformation of the predictor variables allows, first, controlling for collinearity issues (see below) and second, improving the interpretability and comparability of the estimated regression coefficients, especially when interactions are present (Schielzeth 2010).

Accepted Ar ticle
2) We esti ated the coefficients of the selection gradients (the β i , β ij and j ), by fitting the selected best model against the functional traits scaled to mean zero and unit variance.
3) Finally, quadratic coefficients ( i ) were multiplied by two, because the multiple regression procedure gives an estimates equal to 1/2 i (Stinchombe et al. 2008).
Collinearity resulting from correlations between predictor variables is a potential source of problems in multivariate linear regressions and in particular this is expected to affect statistical significance of correlated variables by increasing type II errors (Schielzeth 2010).In our case, with the exception of Dbh and Comp20 (corr =0.7) and LMA and LWC (corr=0.45), the correlations were moderate.Moreover, the variance inflation factors were all below 10 (Supplementary material Appendix1 Table A5, A6, A7, A8).Thus, we expect to keep a relatively high statistical power for parameter significance testing.

Association among functional traits
The first four axes of the PCA explained 67.8% of the total phenotypic variance, indicating a phenotypic structure with non-random association of the eight functional traits (PC1: 22.8%, PC2: 18.8% PC3: 14.1%; PC4: 12.1%).Along the first axis, individuals with early leaf unfolding (i.e.high PSS) and high LWC were opposed to individuals showing high LMA and high d13C (Fig. 2, Supplementary material Appendix1: Tables A2, A3, Appendix2: Fig. A1).Along the second axis, trees with high HV were opposed to those having high LA.The variation in WDens was mainly represented on the fourth axis and orthogonal to N mass in the plan defined by PC3 and PC4.Finally, neither tree height nor competition index was tightly related to any of the four axes (Fig. 2).Based on Table 2, PC1 opposed traits values associated to water uptake (negative PC values) versus traits values associated to water economy (positive PC values).This was also the case of PC2 with negative values associated to water uptake and positive values being associated to water economy.

Effect of competition
The functional traits d13C, PSS, HV and WDens were not significantly associated with competition at any interval of distance (Fig. 3).N mass appeared to be loosely correlated with competition at very short distances (r=0.16,P=0.05).LMA and LA were positively correlated with competition.The strongest relationship was found when considering competitors in a radius of 20 m for LMA (r=0.32;P=1.07.10 -4 ) and of 10 m for LA (r=0.29 P=4.4.10 -4 ).LWC was negatively correlated with competition with the lowest correlation at 20 m (r=-0.23;P= 5.10 -2 ).The performance traits RW10M, FO tree , FO branch were negatively affected by competition, while FW mean was positively related to competition (Fig 3).
Considering the PCA axes (Supplementary material Appendix2: Fig. A2), PC1 was positively correlated with competition, with the highest correlation at 20 m (r=+0.33,P<0.001).PC3 was negatively correlated with competition, with the lowest correlation at 20 m (r=-0.22,P=0.007).PC2 and PC4 were not significantly associated with competition at any interval of distance.

Spatial autocorrelation
The functional traits PSS, LA, d13C, HV, WDens and the performance traits FO branch and FW mean displayed no significant patterns of spatial autocorrelation as shown by the non-significant regression coefficient of Moran's I statistics against pairwise distance (Table 1, Supplementary material Appendix2: Fig. A3B).By contrast, the functional traits LMA, N mass , LWC and the performance traits FO tree and RW10M showed a significant decrease of phenotypic similarity with increasing distance between individuals at the 95% confidence level.However, the distance up to which similarity was higher than expected by chance (D Tresh in Table 1) varied markedly among traits, from 10 m for N mass up to 45 m for LMA.
The three first axis of the PCA showed a significant decrease of phenotypic similarity with increasing distance (Supplementary material Appendix1: Table A4).The distance up to which similarity was higher than expected by chance was higher for PC1 (25 m) and PC2 (30 m) than for PC3 (5 m) (Supplementary material Appendix2: Fig. A3C).

Selection gradients for growth performance
Phenotypic selection analyses at PC-level revealed no significant main directional or quadratic effect of any of the four PCs on radial increment RW10M (Table 3).Moreover, only one interaction was found involving PC2 and PC4 which showed a significantly negative interaction on RW10M, meaning that positive (respectively negative) values on PC2 together with negative (respectively positive) values on PC4 was associated to high RW10M ( PC2-PC4 =-0.08,P=0.02).At trait level, the multivariate model including quadratic and correlational selection gradients (Table 4) revealed that all traits, but LMA, were associated to RW10M either directly or in interaction with the other traits.First, we found negative quadratic selection gradients, suggesting optimal intermediate values maximizing RW10M three traits (LA, N mass , d13C) that apped to independent a es of (co)variation (γ LA =-0.14, P l .5, γ Nmass =-0.19,P .9; γ δ 3 =-0.16,P=0.03).Second, we found significant correlational selection for eleven pairs of traits mostly mapping on different PCs.The observed interactions validated the hypothesis of a synergic effect of traits on performance: for instance PSS and D13C interacted negatively ( PSS-d13C =-0.1, P=0.05) or LWC and WDens interacted negatively ( LWC-WDens =-0.29,P=0.009).However, we also found non synergic effects and for example WDens interacted positively with PSS ( PSS-WDens =0.23, P=0.003).In particular, PC2 traits interacted in a non synergic way with PC1 traits (ex: PSS-LA =-0.11,P=0.04, PSS-HV =0.12, P=0.055).Fitness surfaces were strongly hump-shaped when quadratic selection was significant (Fig. 4A), or were more saddle-shaped when only correlational selection was present (Fig. 4B).Finally, we found a highly significant effect of both dbh and competition on radial incre ent RW M (β dbh =0.77, P< ., β Comp20 =-0.16,P=0.006) (Supplementary material Appendix1, Table A5).Detailed results of phenotypic selection analyses on functional traits can be found in Supplementary material (Appendix 3).
At trait-level, the complete multivariate model revealed that the selection surface for all three reproductive traits was strongly departing from simple linear relationships, with several significant quadratic and correlational parameters estimated (Table 4).First, we systematically found an optimum for a trait belonging to PC1; for PSS in particular, optimums were detected that maximized all three reproductive traits (FO tree : γ PSS =-0.16,P=0.03; FO branch : γ PSS =-0.30,P=0.01; FW mean : γ PSS =-0.24,P=0.05).Besides, optimums were also found for various traits from other PCs and in particular for PC2 traits such as LA for FO tree (γ LA =-0.14, P=0.04) and FO branch (γ LA =-0.34,P=0.002) or HV for FW mean (γ HV =-0.18,P=0.03).As for growth, we observed correlational selection involving all PCs, however we observed more interactions between traits of the same PC and as well as between PC1 and PC2 traits.The interactions involving traits of the same PC tended to have synergic effects (ex: FO tree : LA-HV =-0.09,P=0.045, FO branch : PSS-d13C =-0.22,P=0.003, but see delta13C and LWC), while the opposite was often observed for
For reproductive traits, directional, quadratic and correlational coefficients were of the same order of magnitude, the corresponding fitness surface was slightly humped (e.g.Fig 4C , 4D, 4E).In case of strong directional selection, the fitness surface showed an area of high fitness (Fig 4E).
Finally, competition at 20 m significantly decreased FO tree (β Comp20 =-0.16,P<0.001), while we found no significant effect of dbh on FO tree , FO branch , or FW mean (Supplementary material Appendix1, Tables, A6, A7, A8).Detailed results of phenotypic selection analyses on functional traits can be found in Supplementary material (Appendix 3).

Multivariate phenotypic structure of functional traits
The multivariate analysis shows that the eight studied functional traits are structured in four independent axes that overall explain 67.8% of the phenotypic variation.These observed combinations can be due to a plastic response to environmental heterogeneity such as in soil water content, soil nitrogen content and access to light (Wright and Sutton-Grier 2012) and/or reflect within-stand genetic variation (Donovan et al. 2011).Disentangling the genetic and plastic responses could not be achieved here and was beyond the scope of that paper.Nevertheless, since adaptation is directly related to the phenotype, a reasonable hypothesis is that the observed pattern of phenotypic (co)variation can bring useful insights to adaptation.
The first axis highlighted that individuals with early leaf unfolding had leaf traits allowing water uptake i.e. low water use efficiency, low leaf mass per area and high water content, which is expected to maximise the carbon acquisition through longer effective growing season, jointly with a maximal use of early spring rainfalls before summer drought (Davi et al. 2006).Conversely, individuals with late leaf unfolding also had leaf traits limiting water stress i.e. high water use efficiency, high leaf mass per area and low water content, which is expected to optimize the use of available water resources and would allow copying with the strong summer droughts.Micro-site conditions such as soil spatial heterogeneity (Nourtier et al. 2014) could contribute to the variation in water availability, with trees on deep soils coming with earliest leaf unfolding and higher LWC and trees on shallow soils coming with higher LMA and WUE.Another hypothesis is that phenology itself alters the leaf morphological traits.On the other hand, in a companion study (Bontemps et al. 2016), we found a significant genetic basis for budburst phenology and d13C (less evidence for LMA) which suggests that the presence of genetic correlations cannot be excluded, and that the observed pattern of resource-conservative traits vs resource-exploitative traits can be the consequence of selection.
The second axis opposes trees with small leaves and high branch diameter vs trees with large leaves and small branch diameter.We could also interpret this axis as a water economy/water uptake spectrum, relying on morphological and architectural adaptations instead of phenological/physiological adaptations as for PC1.Leaves with high LA and branches with low HV increase the carbon uptake when water is available, while high HV decreases the risk of embolism during drought.Interestingly, these traits measured at branch scale do not map on the same PCA axis than WDens (correlated with PC4), which corresponds to another organization level (trunk).The variation in WDens is expected to be related to the strong within-population variation in resistance to stem cavitation observed in F. sylvatica (Worteman et al. 2011), even though the relationship between resistance to cavitation and WDens at intra-population level remains to be proven in this species (as in Ruiz-Diaz-Britez et al. 2005 for Pseudostuga menziesii).Finally, N mass behaves quite independently from the other traits and from LMA in particular, in opposition with the observations made by interspecific studies (Wright et al. 2004) but in agreement with other intraspecific studies (Laforest-Lapointe et al. 2014).

Phenotypic selection analyses of functional traits and performance traits
We expected selection to be strong in this site located at the dry limit of the local beech range.Our results supported this hypothesis, as we showed that all studied functional traits, which also involve several uncorrelated traits, were related to the performance traits.Note however that in such in-situ experiment, the correlation between performance traits and the actual fitness is out of reach and part of the variation in performance traits may be neutral.In long-lived organisms, a major issue is the persistence of adaptive growth variation (expected to relate to survival fitness) at the adult stage and whether reproductive performance would provide a more robust evidence of selection.Another issue is the variation of performance across time at adult stage.In this study, we used growth averaged over ten years, assuming that adaptive differences would arise from long-term cumulative effects of climate.However, separating dry and wet years might provide important insights on the role of each functional trait.Regarding seed production, Fagus sylvatica is a mast seeding species and thus, we aimed to capture much of individual fecundity variation within a single mast year.However, we cannot exclude that compensation in non-mast years might occur.
Using various models of selection gradient, we showed that the water uptake syndrome tended to be associated to better performances and especially to higher reproductive outputs.However, we unexpectedly found that for PC2 the water economy syndrome improved performance.The differences among PCs might result from their variation at either the leaf (PC1) or the branch (PC2) scale.On the other hand, our finding that quadratic selection coefficients allowed a better fit also suggested that antagonistic selection was acting on the studied traits.In particular, significant stabilizing selection gradients suggested the existence of optimal trait values for several traits.This was for example the case of leaf unfolding phenology, consistently with the expected trade-off between the increase of the length of the growing season (favouring fruit maturation) and the minimisation of the impact of late frosts.These results agree with our main hypothesis of antagonistic selective pressures resulting from the maximisation of carbon acquisition versus the capacity to cope with strong summer droughts (Lo Gullo and Salleo 1988).

Accepted Ar ticle
One major issue in this study relates to the possibility that small-scale environmental variation may result in the joint phenotypic response of both functional and performance traits to environmental variation, in particular for the first and third axes, making it difficult to establish causal links between functional and performance traits variation and firmly conclude on selection rather than simple environmental covariance (Rausher 1992).In our case, we observed better growth and reproduction associated to traits favouring resource acquisition in agreement with a positive effect of resource availability on tree performance.On the other hand, we were also able to model more complex selective surfaces showing that trait values associated to a priori less favourable environments can also promote tree performance.
Finally, our results suggested that several adaptive strategies may coexist within the studied natural population.The existence of interactions effects among functional traits, showing that the relationship between a trait and performance depended of another trait value, suggests that the combinations themselves affect tree performance.For example, we found that individuals with early leaf unfolding had the highest reproductive outputs provided that they also had strong leaf water content.was not the case for late bursting individuals.Such interaction suggested a synergic effect of the traits.However, we also found many interactions, suggesting that synergy might not be systematic, especially when interactions involved traits varying at different scales i.e. between traits from PC1 and PC2.
Overall, our results suggest that various combinations of the water uptake strategy and the water economy strategy co-exist within our population.This means that the impact of each trait on performance highly depends on the integrated phenotype and this, together with the previous findings, strongly supports the importance to assess the adaptive value of coordinated suites of traits (Laughlin and Messier 2015).

Table Legends
Table 1: Variation of the functional traits (Func), traits used as covariates (Cov) and performance traits (Perf).N is the number of trees successfully phenotyped (among the 149 targeted).Results of spatial auto-correlation analyses are given by D tresh (the distance up to which similarity is higher than expected by chance), and Slope (the slope of the regression between trait similarity and distance).

Figure 2 :
Figure 2: Multivariate phenotypic correlation structure characterised by a principal component analysis on eight functional traits (PSS: leaf unfolding phenology; LMA: leaf mass per area; LA: leaf Area; HV: Huber Value; WDens: wood density; LWC: leaf water content, N mass : leaf nitrogen content; d13C: leaf carbon isotopic discrimination).Tree height (Height) and competition intensity from neighbours located at a distance up to 20m from the target tree (Comp20) were included as supplementary variables.The first and second axes (a) and the third and four axes (b) are plotted.Accepted Ar ticle

Figure 3 :
Figure 3: Correlation coefficients (r) between the Martin-Ek competition index at a given distance and leaf functional trait values.See table 1 for the codes of the functional traits.Competition intensity was estimated considering the potential competitors within a circle with radius ranging from 1 m to 20 m and a 1 m step.Accepted Ar ticle

Figure 4 :
Figure4: Fitness surfaces resulting from the combination of directional, correlational and quadratic selection gradient for selected pairs of traits and different proxies of fitness: A, B: RW10M; C,D: FO tree E: FO branc F: FW mean .All the fitness surfaces shown correspond to significant interaction effects between traits, but in A-C-D-E-F quadratic terms were also included in the model, while in B, there were no other effects.