Uniform Selection as a Primary Force Reducing Population Genetic Differentiation of Cavitation Resistance across a Species Range

Background Cavitation resistance to water stress-induced embolism determines plant survival during drought. This adaptive trait has been described as highly variable in a wide range of tree species, but little is known about the extent of genetic and phenotypic variability within species. This information is essential to our understanding of the evolutionary forces that have shaped this trait, and for evaluation of its inclusion in breeding programs. Methodology We assessed cavitation resistance (P 50), growth and carbon isotope composition in six Pinus pinaster populations in a provenance and progeny trial. We estimated the heritability of cavitation resistance and compared the distribution of neutral markers (F ST) and quantitative genetic differentiation (Q ST), for retrospective identification of the evolutionary forces acting on these traits. Results/Discussion In contrast to growth and carbon isotope composition, no population differentiation was found for cavitation resistance. Heritability was higher than for the other traits, with a low additive genetic variance (h2 ns = 0.43±0.18, CVA = 4.4%). Q ST was significantly lower than F ST, indicating uniform selection for P 50, rather than genetic drift. Putative mechanisms underlying QST


Introduction
The climatic niches of forest tree species are moving faster than the maximum rate of migration, measured by palynology or gene flow analysis, as a direct consequence of the current increase in temperatures due to global warming [1,2]. Forest tree populations are thus facing new selection pressures and are unable to track their bioclimatic envelope [3] over the time scale at which these changes are occurring. The local adaptation and survival of tree populations in a rapidly changing environment with warmer temperatures and more frequent water shortage is a major concern in efforts to ensure the sustainability of forest ecosystem services. In addition to this trend, climate experts are predicting more extreme climatic events, such as periods of severe drought [4], which will increase mortality rates [5,6,7]. These effects on tree mortality highlight the way in which the impact of climate change may depend on the changes associated with extreme events rather than trends [8]. In this context, there is a need to investigate relevant drought tolerance-related traits, to quantify both genetic variation and phenotypic plasticity, which together define the capacity of tree populations to adapt.
From a physiological point of view cavitation resistance is an important trait to estimate drought tolerance. Indeed, dysfunctions of the vascular system of the tree, such as xylem embolism due to cavitation events, is likely to be a key factor governing the mortality of these long-lived organisms [9].When a cavity is formed in the xylem sap under tension (negative pressure), it may spread in the vascular system through intervessel or intertracheid pits, thus compromising the capacity of the plant to transport water [10].
A direct causal link between survival (fitness) and cavitation resistance during extreme drought has been highlighted, based on two lines of evidence suggesting that cavitation resistance is an important adaptive trait. Firstly, assessments of the correlation between cavitation resistance and lethal water potential [11,12] demonstrated a highly significant linear relationship (r 2 = 0.9) between these two traits. Secondly, global surveys of cavitation resistance in woody species showed that xeric species are more resistant to embolism than hydric species [13,14,15]. These interspecific studies, with adaptive inferences concerning cavitation resistance being rendered robust by the incorporation of phylogenetic information [15,16,17], concluded that cavitation resistance-related traits are under natural selection [15,18]. To validate this evolutionary pattern, a population-level perspective is appropriate, because variation observed across species cannot be assumed to reflect patterns within species.
At the intraspecific level, cavitation resistance can be analyzed by provenance or progeny trials. The few studies carried out to date (reviewed in Table S1) have included only small numbers of individuals (,9) and populations (,5), and it has therefore not been possible to estimate environmental and genetic effects on phenotypic variation accurately. We therefore still know little about the genetic determinism and micro-evolutionary pattern of this hydraulic trait, but such information is absolutely necessary if this trait is to be incorporated into breeding programs and for a more fundamental understanding of the evolutionary basis of tolerance to severe drought.
The aim of this study was to provide the first estimates of heritability, additive variance and population differentiation for cavitation resistance-related traits. We carried out a case study of maritime pine (Pinus pinaster Ait.), a species with a fragmented distribution in the western part of the Mediterranean region. The scattered distribution of this species may have prevented or limited gene flow between different groups of populations, promoting high levels of genetic divergence between ecotypes due to genetic drift [19] and/or natural selection (Quezel and Barbero 1998 in [20]). Here, we took advantage of a new technology (high-throughput phenotyping platform for cavitation resistance) to screen for the first time a large number of genotypes from six ecotypes to test the hypothesis that Pinus pinaster populations have been subjected to diversifying selection for cavitation resistance. More specifically, this experiment aimed to address the following questions: what is the level intraspecific variation and heritability for cavitation resistance? Can we separate the relative roles of drift and selection in population differentiation for this trait?

Provenance trial and climatic data
We carried out a provenance-progeny trial, in which young trees (six-year-old plants) were planted in December 2003 at the INRA forestry station in the Aquitaine Region (44u44'N, 00u46'W). The mean annual temperature at this site is 13.2uC and mean annual rainfall is 836 mm . The soil is a sandy podzol with a water table rising to about 0.5 m below the surface in winter and descending to a depth of 2 m in late summer. Seedlings were grown in the nursery from open-pollinated seeds collected from 24 natural populations (or ecotypes) in France, Spain, Morocco and Tunisia, to cover the fragmented distribution of Pinus pinaster. Each population was represented by 20 to 30 halfsib families. The trial was arranged in a randomized block design (15 blocks) with single-tree plots. Each block contained at least one tree from each half-sib family. There were 600 seedlings per block, giving 9,000 seedlings for the entire experiment.

Choice of populations
The assessment of cavitation resistance is a intensively timeconsuming process [21]. We therefore designed a procedure for the selection of a subset of populations representing all the variability of climatic envelope of maritime pine. For a total of 754 grid points covering the entire natural range of the species [22], we first extracted climatic data from the CRU CL 2.0 109 global dataset for the period 1961-1990 [23,24,25]. These data included monthly average precipitation, mean, minimum and maximum temperature, diurnal temperature range, water vapor pressure, cloud cover, wet day frequency, ground frost frequency, radiation, wind, Martonne index, Turc's potential evapotranspiration and soil water deficit. We also derived the vapor pressure deficit from these parameters. Principal component analysis (PCA) was then used to reduce the number of dimensions for the whole set of climate variables ( Figure 1). The data were centered and scaled before PCA. The 14 populations were finally placed on the main plane of the PCA (accounting for 76% of the variation, Table S2) and six of these populations (Table 1) were selected as a representative subset of the climatic envelope covered by Pinus pinaster species. In each population, eight families (5 half-sibs/family/block) were randomly selected for further analysis.

Sample preparation for the assessment of cavitation resistance
We collected branches, according to the sampling procedure described below, during winter 2009, before 10 am, to avoid native embolism. Needle water potential was lower than -1 MPa, far from the minimum needle water potential in summer (22 MPa) of Pinus pinaster [26]. The branch sample corresponded to the 2007 and 2008 growth units on the 2007 whorl when possible, in order to measure on the same number of rings in each sample. Sampled branches were fully exposed to the sun, longer than 40 cm and with a diameter of 0.3 to 1 cm (,4 years of age). The current needles were removed and the branches were wrapped in wet paper towels and bagged upon collection to prevent dehydration. In the lab, samples were cut under water just before measurement to obtain 0.28 m segments (i.e., much longer than the longest tracheid). Bark was removed from all segments before measurements.

Assessment of cavitation resistance
Cavitation resistance was measured on 240 genotypes (6 populations * 8 families * 5 offsprings), with the Cavitron technique [21,27,28]. Centrifugal force was used to establish negative pressure in the xylem and to provoke water stress-induced cavitation, using a custom-built honeycomb rotor (Precis 2000, Bordeaux, France) mounted on a high-speed centrifuge (HS18, MSE Scientific, London, UK). This technique enables measurement of the hydraulic conductance of a branch under negative pressure. Xylem pressure (P x ) was first set to a reference pressure (20.5 MPa) and maximal conductance (k max ) was determined by measuring the flux of a reference ionic solution (10 mmol dm 23 KCl and dm 23 mmol dm 23 CaCl 2 in deionized water) through the sample. The centrifugation speed was then set to a higher value for 3 min to expose the sample at a more negative pressure. Conductance (k i ) were measured 4 times for each step, and the average was used to compute percent of loss of xylem of conductance (PLC in %) following PLC = 100 (12k i /k max ). The procedure was repeated for at least eight pressure step with a 20.5 MPa step increment until PLC reached at least 90%. The percent loss of xylem conductance as a function of xylem pressure (MPa) represents the sample's vulnerability curve ( Figure 2). Rotor velocity was monitored with a 10 rpm resolution electronic tachymeter (A2108-LSR232, Compact Inst, Bolton, UK) and xylem pressure was adjusted to about 60.02 MPa. We used Cavi_soft software (version1.5, University of Bordeaux) for conductance measurements and computation of all vulnerability curves (VCs). The 10,800 measurements of conductance were performed at the new high-throughput phenotyping platform for hydraulic traits (CavitPlace, University of Bordeaux, Talence, France).
Based on a sensitivity analysis (data not shown) and graphical checking (Figure 2), we retained the reparameterized sigmoid function fitted to the conductance data (k i ) (see [29] and [30] for an exhaustive review) rather than the Weibull model, for determination of the pressure at which the sample lost 50% of its conductance (P 50 ) and the slope of the curve at P 50 (S 50 ). The contour plot represents the presence's probability (kernel density estimate) of Pinus pinaster population (small black dot) within the bioclimatic envelope representing by PC1 and PC2 axes, accounted for 54% and 21% of the variance, respectively. The studied populations and provenance test are indicated by black circles. PCA was performed with the variables indicated in the methods section. See Table S2, for additional information about the relative contribution of climatic variables to the axes. Lower panel: projection of 14 climatic variables on the subspace spanned by the first two eigenvectors (correlation circle). doi:10.1371/journal.pone.0023476.g001 where k max is the highest conductance measured for each sample (equivalent to k sat in the original Ogle's model), k i is the mean conductance at a given pressure concerned, X is the percentage loss of conductance of interest (in %), P X (in MPa) is the pressure inducing X% loss of conductance and S X (in MPa.% 21 ) is the slope of the tangent at the P X abscissa point on the curve. Analysis has been performed for pressures and slopes corresponding to X = 12, 50 and 88 % loss of conductivity (P 12 , P 50 , P 88 , S 12 , S 50 and S 88 respectively).

Carbon isotope ratio and growth measurement
Carbon isotope ratio (d 13 C in %o) was obtained as previously described [32,33]. Needles of the growth unit used for cavitation analysis were harvested and 20 needles were sampled at random. The needles were dried and ground to a powder and 3 mg sample was analyzed with an isotope ratio mass spectrometer (FISONS Isochrom, Manchester, UK) at INRA facility in Reims (France). Total height was measured at the ages of two (2004) and three (2005) years, on the same six populations and eight families, for all 15 blocks. The annual increase in height (D h ) was calculated as the difference between these two measurements (in 2004 and 2005).

Quantitative genetic analysis
Genetic analysis was conducted with the following mixed model: where y is the vector of observation for a trait, b is the vector (number of block) of fixed block effects, pop is the vector (number of populations) of random population effects, f is the vector (number of mother trees) of the random genetic effects of mother tree within the population, e is the vector (number of individuals) of residuals, X is called the design matrix, Z 1 and Z 2 are the incidence matrices linking the observations to the effects. A variance was fitted for each random effect: s 2 pop is the genetic variance between populations, s 2 f (pop) is the genetic variance between mother trees nested within a population and s 2 e is the residual variance.
Variance or covariance components were estimated by the restricted maximum likelihood (REML) method, assuming a normal distribution of the random effects. The significance of variance components were tested using log-likelihood ratio tests. We included population as a random effect to draw inference at species levels [34] and to obtain an unbiased estimate of heritability and genetic population differentiation [35]. The normality, identity and independency of residuals of each trait were graphically checked by plotting studentized marginal and conditional residuals (available on request), which confirmed that the data match with the assumption of mixed model. We estimated narrow-sense heritability as follows: where s 2 A is the within-population additive variance. In our study, s 2 A was estimated by s 2 A = 4s 2 f (pop) as trees from the same family were presumed to be half-sibs (open-pollinated seeds). We did not include the population effect in the heritability calculation, because natural selection appeared to occur within each population [36]. The standard deviation of heritability was calculated with the equations of delta method (see Appendix in [37]).
Variance components were standardized by the trait mean [38] as follows, CV X = 100!(Variance)/Mean X where X is the trait considered, and CV is the coefficient of variation. Each variance component is expressed with a CV (CV A : additive coefficient of variation; CV BP (V BP = s 2 pop ): coefficient of variation between populations; CV P : phenotypic coefficient of variation; CV R : residual coefficient of variation). The variance of each component was extracted from the asymptotic covariance matrix. The significance of mean population difference was estimated using the same model (Eq. 2) with a proc GLM with a Student-Neuman-Keuls post hoc test.

Correlation between traits
To facilitate interpretation of correlation, negative value of P 50 and d 13 C were converted from negatives to positives. Genetic correlations between traits were evaluated by calculating Pearson's coefficient on the family Best Linear Unbiased Predictor estimation (BLUP, for additional information see [37] p745). BLUP estimation ensures that data are corrected for block effect. We will refer to these correlations as genetic correlations. For phenotypic correlation, all Pearson correlations were computed over the BLUP family plus BLUP population and the grand-mean. Estimation of population differentiation The estimate of phenotypic differentiation between populations, Q ST [39], was calculated as Putatively neutral nuclear microsatellites (nuSSRs) were used to account for genetic differentiation (F ST ) caused by demographic and other processes not related to selection (e.g., genetic drift resulting from geographic isolation or population expansion). Eight polymorphic nuSSRs were selected from those previously developed by [40] (NZPR413, NZPR1078, ctg64), [41,42] (ctg275, FRPP91, FRPP94, ITPH4516) and [43] (A6F03). The markers were selected to be evenly distributed over the various linkage groups of the maritime pine genetic linkage map, with, at least, 4 alleles and multiplexing capacity. Genotyping was performed on genomic DNA isolated from the needles of 20 to 30 individuals from each of the six selected populations, as previously described [44,45]. We used F ST (which is estimated from the allelic frequency) rather than R ST (which also takes into account allele size) because genetic drift affects allele frequency but not mutation rate. F ST values for each locus were estimated with Genepop [46] using the framework developed by [47] adapted for SSR data [48]. For the comparison of Q ST and F ST , to disentangle the effects of genetic drift from those of selection, we develop a new test to avoid previously reported limitations [49,50,51,52] and allow to test Q ST .F ST and Q ST ,F ST . We explicitly derive the F ST and Q ST distribution using in both case a parametric bootstrap as follows.
F ST * is a parametric bootstrap replicate of F ST . First, nuSSR loci were randomly resampled with replacement, to estimate the sampling variance of F ST . Each of this F ST replicate (F boot ST ) value was then multiplied by a random number drawn from the Lewontin-Krakauer distribution, a chi-squared distribution with a number of degree of freedom equal to the number of loci minus 1, divided by degree of freedom equal to the number of loci minus 1. This distribution has been shown to take into account most of the deviation from the neutral model due to demographic history [51,53]. We will refer to this distribution of F ST * as the ''drift distribution''.
We estimated the sampling variance of Q ST (Eq 4), by simulating the distribution of each variance component (s 2 pop ,s 2 A ) with a parametric bootstrap [50], using the Satterthwaite approximation [54]. This distribution is highly conservative and takes into account the deviation from homogeneity of variance [54].
dfe i is the effective degree of freedom for the variance component of i factor. df i is the observed degree of freedom. s i is observed variance component due to i factor. s i{1 is observed variance component due to i-1 factor (nested or residuals factor). n is the total size of the sample. j is the number of level of factor i. We calculate a Q ST * for each replicate from s 2Ã pop and s 2Ã A (Eq. 6). We will refer to this distribution of Q ST * as the ''phenotypic distribution'', although Q ST is a standardized measurement of additive genetic variance between populations.
Finally, we compared the F ST * and Q ST * distributions, using nonparametric and free distribution two-sample test for equality of the 2.5 and the 97.5 quantile (see [55] for theoretical proof) with a Bonferroni correction for multiple comparisons. We also applied a studentized bootstrap, which gave similar results but required more computation time [56,57,58]. All the analyses were performed with SAS version 9.2. Codes are available on request.

Between-population variation
For each population, vulnerability curves showed similar sigmoid shape with the air-entry (P 12 ) around 23.2560.006 MPa (see Figure 2). Linear curves were discarded from the analysis [59]. The between-population effect (V BP ) was significant for d 13 C and D h but not for P 50 (Table 2). Similarly, no difference was found for the other cavitation resistance-related traits (S 12 , S 88 , S 50 , P 12 , P 88 , data not shown). Cavitation resistance-related traits had much lower coefficients of variation than D h (Table 2). This was particularly true for the between-population coefficient of variation (CV BP = 1% and 18% for P 50 and D h , respectively). It should be noted that CVs for d 13 C are not comparable with those of other traits, because they are estimated relative to a standard [60] and are therefore independent of scale change but not of origin. The fixed block effect was significant for all the traits studied, indicating that some of the environmental variation was taken into account by the experimental design.
The populations from the wettest areas (Mimizan and San Cipriano) had the highest D h values (Figure 3a), whereas Tamrabta population (from Morocco) presented the lowest value. Iberian populations from very different climatic areas (Coca, Bayubas, Oria) had intermediate values, with no detectable trend as a function of environmental aridity.
No significant difference between populations was detected for P 50 (Figure 3b), although Tamrabta surpassed the other populations and was the most cavitation-resistant population. Tamrabta also presented the lowest d 13 C value (Figure 3c), demonstrating a significantly lower water-use efficiency than the other populations, all other populations presented similar d 13 C values.

Within-population variation
Heritabilities and normalized measurements of trait dispersion (i.e. CVs) were estimated to evaluate the within-population additive variance, evolvability (through the analysis of CV A ) and micro-environmental sensitivity (through the analysis of CV R ) ( Table 2). Narrow-sense heritability (h 2 ns ) for P 50 was higher (0.4460.18) than those estimated for d 13 C (0.2160.10) and D h (0.3560.06), showing that cavitation resistance was genetically controlled, although the standard error was high, probably due to the small number of progenies per mother tree analyzed. The CVs of P 50 and D h presented contrasting patterns, with a lower coefficient of additive variation for P 50 (CV A = 4.4%) than for D h (CV A = 16.2%), suggesting limited evolvability of P 50 .

Correlation between traits
We found a significant positive phenotypic correlation between absolute value of P 50 and d 13 C at the phenotypic level (r = 0.30, P = 0.035 see Figure S1), indicating that the more cavitationresistant genotypes tended to be less water-use efficient. However, this relationship was not significant at the genetic level (r = 0.14, P = 0.320). No relationship between P 50 and D h was found at either phenotypic or genetic level. A significant negative phenotypic correlation between D h and d 13 C was detected (r = 20.68, P,0.0001). This correlation was barely significant at the genetic level (r = 20.29, P = 0.053).

Discussion
We reliably estimated for the first time the genetic variability of cavitation resistance, a functional trait that allows plants to survive under severe drought. We also provided evidence of natural selection acting on this trait. These results were based on the greatest number of genotypes ever measured to date in an experimental design (Table S1 and supplementary references). Despite the high level of variation of cavitation resistance between species [15], we detected no significant differences between maritime pine populations from a wide range of environments. Moreover, the between-population variability of cavitation resistance was significantly lower than would be expected under a hypothesis of genetic drift alone (Q ST -F ST comparison). We can therefore reject the hypothesis of diversifying selection. We suggest instead that uniform selection has shaped the phenotypic variability of this trait. Uniform selection could be seen as a stabilizing selection acting within each population with the same selection optimum in each population despite the steep climatic gradient [62,63]. Conversely, growth and water-use efficiency displayed different patterns and were found to be subject to strong diversifying selection and genetic drift, respectively. Quantitative genetics analysis also showed that cavitation resistance presented a significant heritability, higher than that estimated for growth and water-use efficiency. This is the first evidence of uniform selection in woody plants and the underlying mechanisms are discussed below from a micro-evolutionary point of view.

Intra-vs. interspecific variability of cavitation resistance
Despite the steepest climatic gradient (precipitation ranging from 400 to 1,600 mm in the sampled populations) and strong phylogeographic structure between the six studied populations [22,64]. Very low between-population variance for cavitation resistance were found (CV BP = 1%). The few studies published to date (reviewed in Table S1) tended to skim over the issue of intraspecific variation of cavitation resistance in provenance or progeny trials and reported little or no difference between populations [65,66]. As these studies were not designed to assess the genetic component of phenotypic variation, further investigations are required, to generalize our finding to other species. In addition, phenotypic variation for cavitation resistance was low (CV P = 6.6%), but consistent with the range reported for wood properties of maritime pine, such as mean ring density, lignin content and fiber morphology ( [67,68] Lamy, unpunblished data). However, we are lacking information about the intraspecific variation of pit pair anatomical traits that are known to be implicated in cavitation resistance [69,16]. The low within-species variability for cavitation resistance is remarkable (P 50 = 23.9360.04 MPa, estimated over the whole dataset), given that substantial variability has been described between species. For instance, Delzon et al (2010) showed that cavitation resistance ranged from 23 to 212 MPa in a sample of 40 coniferous species. This variability was interpreted as the effect of natural selection rather than phylogenetic legacy [15].
In contrast, the population differentiation observed for growth and water-use efficiency (WUE) was significant and consistent with previous results for this species [70]. The Moroccan population had the lowest WUE, consistent with previous findings based on both gas exchange measurements and carbon isotope discrimination [71,72,73,74]. In a provenance trial carried out in south-western France, this Moroccan population displayed lower stomatal sensitivity to water stress (delayed stomatal closure), leading to greater water loss throughout the summer period and a lower WUE (as reflected by carbon isotope composition). Genotype 6 environment interaction could potentially alter differences between populations [75,76]. Our results therefore require confirmation in provenance trials carried out in drier climates.

Relationships between traits
The weak but significant positive correlation found between absolute value of cavitation resistance and water-use efficiency (carbon isotope composition) suggested that drought-tolerant genotypes had lower water-use efficiency. In dry environments, genotypes that allocate more carbon to the construction of cavitation-resistant wood in order to avoid runaway embolism might be able to maintain higher stomatal conductance and hydraulic conductance at low leaf water potential, resulting in a decrease in water-use efficiency. Our results are consistent with previous findings [77] of a strong and positive relationship between these two traits in two cedar species. However, little or no correlation has generally been reported [78,79,80]. The Table 2. Variance components (V P , V BP , V A , V R ), narrow-sense heritability (h 2 ns ), coefficient of variation (CV P , CV A , CV BP , CV R ) and population differentiation (Q ST ) for all studied maritime pine populations. h 2 ns is the narrow-sense heritability and SE is the standard error of heritability, V P is the phenotypic genetic variance, V A is the additive genetic variance, V BP is the between-population variance, V R is the residual variance. CV A is the variation coefficient of additive variance after adjustment for the block effect. CV P is the variation coefficient of phenotypic variance after adjustment for the block effect. CV R is the residual coefficient of variation. CV BP is between-population coefficients of variation. Q ST is the genetic quantitative variation between populations (Spitze, 1993). The significance of random effects is indicated after each variance estimator: ns P . 0.05, * P,0.05, ** P,0.01, *** P,0.001. a CVs for d 13 C are not comparable with other traits as they are estimated relative to a standard. doi:10.1371/journal.pone.0023476.t002 negative correlation between carbon isotope composition and growth has been reported in previous studies [77,81,82], assuming that growth is a function of carbon assimilation and carbon isotope composition is an index of retrospective gas exchanges.

Evidence of uniform selection for cavitation resistance
The phenotypic distribution of cavitation resistance was significantly lower than the expected distribution under the drift hypothesis. This may be interpreted as a consequence of uniform selection (also called homogenous, spatially homogenizing, con-vergent selection, uniform stabilizing selection or stabilizing selection across population). This inference (Q ST ,F ST ) may result from an underestimation of Q ST variance [52], leading to a false positive result. However, the F ST estimate was more than five times greater that the Q ST estimate. We thus believe that this difference is biologically meaningful and not due to a statistical artifact [52]. The robustness of this result is, also, supported by fourth lines of evidence: (i) the different patterns obtained for growth and WUE in the same experimental design. (ii) Selection procedure of population (see methods) increase the probability to find diversifying selection because we selected extreme populations in term of climatic origin and evolutionary history (different mitotypes and chlorotypes, see [64]), consequently uniform selection could not be interpreted as a sampling bias. (iii) Wood density (measured by X ray on the same data set) showed exactly the same pattern (Lamy, unpunblished). (iv) Willson et al (2008) showed from interspecific data with narrow taxon sampling (limited to the Juniperus genus), that cavitation resistance gave strong phylogenetic conservatism, suggestive of uniform selection for the maintenance of ancestral traits.
For growth, diversifying selection was highlighted by Q ST being greater than F ST [83]. For WUE, we detected no signature of selection, as the phenotypic and drift distributions did not differ significantly, as previous reported by [84] for Quercus suber (P ST / F ST comparison). However, the distribution (Q ST -F ST ) of this trait was shifted to the right (integral probability above 0 . integral probability below 0, Figure 4c right panel), suggesting that diversifying selection was slightly more pronounced than drift.
What are the mechanisms behind ''uniform selection'' for cavitation resistance?
The causal mechanisms underlying uniform selection or leading to evolutionary stasis are not well understood [63,85,86,87,88]. We discuss here only the processes most likely to account for the observed pattern (Q ST ,F ST ).
Weak molecular variation. Lethal or sublethal mutations may limit the variability of the genes they affect, thereby controlling trait variation. In conifers, cavitation resistance is known to be determined by xylem anatomy, including, in particular, the characteristics of intertracheid pits [16,69]. As knowledge about the nucleotide diversity of different functional categories of genes accumulates, it may become possible to test the hypothesis that genes involved in intertracheid pit formation (once these genes have been identified) display lower levels of diversity.
Genetic constraints. If selection acts on a trait that is negatively correlated with another trait (or traits) also under selection, than the decrease of rate of evolution for the first trait is proportional to the strength of the correlation [94]. A multi-trait approach could be used to explore this hypothesis indirectly [95], but could fail if the trait is canalized.
Canalized trait. This hypothesis suggest that cavitation resistance is canalized to buffer the variation of this key hydraulic trait against all kinds of disturbance, being of genetic (mutation, hybridization, recombination) and/or environmental nature [89,90,91]. Emergent properties of molecular networks could buffer molecular variability [92], to maintain phenotypic function, in accordance with the robustness theory [93]. In zoology, dipterian wings shape or centroid size (or mammalian body temperature) are the best known cases of canalized traits [90]. Indeed they reported a similar wing shape between species despite a great variability of climatic niche and a low additive genetic variance between populations for this trait. Except for leaf shape in Arabidopsis thaliana, there is no evidence of canalized trait in plants nowadays. For cavitation resistance, two arguments lead  us to consider canalization as the most likely mechanism: (i) low additive genetic variance between populations (V BP ) and (ii) the similarity of cavitation resistance values among all the Pinus species [69].

Variance component analysis
Variance component analysis for cavitation resistance resulted in the first estimates of the heritability of this trait (h 2 ns = 0.4) and its additive coefficient of variation (CV A = 4.4%). These values suggest that this trait may respond to truncation selection frequently practiced in breeding. However, for a given selection intensity, genetic gain for cavitation resistance would be limited by the low additive variance, although long-term artificial selection experiments have shown that quantitative traits have a non negligible mutational variance [91,96,97], which could supply further additive variance at each generation. However, due to the small number of half-sib families and progenies within each family, which could inflate the value for heritability [37], this estimate should be interpreted with caution. Additional studies with a larger sample size are required for further exploration of the genetic determinism of this hydraulically important trait.
The much higher CV A value (16.2%) for height increment is consistent with previous reports and accounts for the genetic gain achieved for this trait over successive generations in breeding programs [98,99,100]. For d 13 C, a previous study [31] reported a slightly lower heritability (h 2 ns = 0.17 vs. 0.21 here), but with estimation based on a diallel cross of limited size, with a narrow genetic background restricted to 12 elite trees from south-western France.

Future directions
Phenotypic variation is a fundamental prerequisite for evolution because natural selection acts on phenotype. Adaptation and evolution via natural selection requires the presence of genetic variation among individuals in a population upon which natural selection can act. Intra-population genetic variability can thus be seen as the fuel for future adaptation. However, an environmentally induced shift in phenotype is also a major component of the variation we see in nature. Recent studies [80,101] showed a weak but significant phenotypic variability for cavitation resistance. Our results suggest that this between populations variability might be under environmental control rather than genetic determinism. These considerations call for more research (ongoing) aiming at quantifying the in situ phenotypic variability of cavitation resistance and the extent of phenotypic plasticity using provenance trials installed under different edapho-climatic environments. Further studies are also being pursued to dissect the genetic architecture of cavitation resistance to determine the number, map location and effects of Quantitative Trait Loci controlling part of the variation of this trait. Figure S1 Genetic (right panel) and phenotypic (left panel) correlation between traits. For ease of interpretation, we have converted all the negative values to positive values (P 50 , dC 13 ). For the genetic correlation, all Pearson correlations (r) were computed over the best linear unbiased prediction (BLUP) (n = 48 for P 50 , d 13 and n = 151 for Dh). For phenotypic correlation, all Pearson correlations were computed over the BLUP family plus BLUP population and the grand-mean, to ensure that the order of degree of freedom remained the same and the block effects are removed. P 50 , pressure at 50 % loss of conductivity in MPa, Dh the annual increment between 2004 and 2005, in mm, d 13 C is the isotope discrimination for carbon 13 in %.

(TIF)
Table S1 Review of intraspecific studies for cavitation resistance estimated using P 50 or related parameters (as indicated in the table). Npop: number of populations used, Nind: number of individuals per population used to assess cavitation resistance. The table is divided in two parts, the first part corresponds to provenance or progeny trials, and the second to ''in situ'' studies. W is mean wet ground days (days). I is mean Martonne's index (P i /(T a +10)). P i is the mean precipitation (mm.days 21 ). C is percent of cloud cover (%). S is the mean of wind speed (m.s 21 ). V is the water vapor pressure in air (hPa). VPD is the water vapor pressure deficit of air (hPa). T min is the minimum temperature (uC). D DT is the mean diurnal temperature range (uC). T m is mean temperature (uC). R G is mean global radiation (W.m 2 ). T max is the maximum temperature (uC). H is mean soil water deficit (P i -ETP, in mm). ETP is mean Truc's potential evapotranspiration (mm). (DOC)