Plant material and seed collection
We sampled 22 mature trees (see below for details on tree selection) from a natural maritime pine stand located in a continuous forest in Central Spain (Coca, Segovia). No forest management has been applied to this area in the last 40 years except for removal of dead trees. Mature, well-formed cones were collected from each tree for progeny testing. The cones were collected for the same year, to avoid interannual environmental effects, and for each mother tree separately. Thus, each seed lot had an approximate structure of half-sib families (since correlated paternity within maternal progenies is known to be negligible in maritime pine, see [30]). Seeds were extracted from several cones from each mother and stored in dry conditions at 4°C.
Reproductive success and field measurements
Twenty-two mother trees were selected to cover the variation in female reproductive success found in the stand. Previously [26], a study on gene dispersal and female reproductive success based on parentage analysis of naturally-regenerated saplings using nuclear microsatellite (nuSSR) molecular markers had been conducted. This study included 380 mature trees sampled in an area of 3.2 ha and 267 naturally-established saplings, considered the basis for next generation’s population, collected from a central subplot of 30 m in radius. Female parentage LOD scores were retrieved from this study. Briefly, the most-likely parents and parent pairs were detected using log-likelihood ratios (LOD-scores) [31],[32] and population allele frequencies estimated from the whole dataset (n = 581). Only highly reliable parent-offspring matches (LOD-score > 4.4, as shown by simulation) were used (see details in [32],[33]). Following González-Martínez et al. [26], the parent closer to the offspring was considered to be the mother in parent pairs, also because pines present a much more restricted seed than pollen dispersal. LOD scores for each mother-offspring relationship were summed across all offspring for each mother and then divided by the sum of LOD scores across all mothers, producing a relative measure of female reproductive success for each mother tree. Estimates of reproductive success based on fractional paternity assignment are normally less biased than those relying on direct parent-offspring identification [34]. To remove biases due to distance effects, individual reproductive success values were corrected by the probability of a propagule dispersed from the focal female landing within the central subplot where saplings were sampled (assumed for mathematical convenience to be a square of equal area than the actual central subplot). This probability was calculated via numerical integration of a bivariate exponential-power dispersal kernel with parameters α = 68.6 and β = 3, which are the best out of three inverse-modeling dispersal kernels for saplings fitted in [26].
For each of the 22 mother-trees with contrasting values of relative reproductive success (see Additional file 1: Table S1 in Supplementary Information for a summary), we measured different fitness-related traits: size (measured by total height –HT– and diameter –DBH–), diameter increment in the last 10 years (DBH10), number of cones, seed weight (based on a 100-seed sample), germination rate (after 35 days), and number of empty seeds (in a 100-seed sample).
Experimental design for progeny testing
To estimate breeding values, two different experiments were conducted using progenies from mother trees with contrasting effective reproductive success. The first one (Experiment 1) was a growth-chamber experiment under environment-controlled conditions, the second (Experiment 2), an outdoor experiment under semi-natural conditions. Individual identification of the seeds and seedlings was maintained for the duration of the experiments.
Experiment 1: Environment-controlled conditions
A sample of 100 seeds from each of 18 mother trees (see Additional file 1: Table S1 in Supplementary Information) were weighed and germinated under controlled conditions (20°C, 12 h light) using a sand-substrate. After germination, 24 seedlings per family were transplanted to individual 400°Cc pots, with a sand-peat (2:1) substrate, and were moved to a growth chamber under controlled conditions. Prior to transplanting, each pot, including the substrate, was weighed (dry weight and full capacity weight) to determine water content at full capacity. The plants were arranged in four blocks (two for each watering regime, see below) and a row-column design (within blocks), with six experimental units (one-seedling plot) for each genetic entry, hence totaling 24 seedlings per family. The conditions of the experiment were set to hasten ontogenetic changes [27], aiming at maximizing developmental differences in a short cultivation time. This, in addition, helped to control root restriction effects due to container to a minimum [35]. The experimental phase was divided in three (acclimation, growing –including the watering treatment period–, and hardening), depending on temperature and photoperiod (see Additional file 1: Table S2 in Supplementary Information). The first watering regime consisted on keeping plants at full-field capacity (W1, 45% in volume), and the second on reducing water availability to 20% of full-field capacity (W2, 10% in volume). The two watering regimes were controlled by weighing two pots per family and block every second day. In total, the experiment lasted for 189 days (84 days under watering treatments).
Experiment 2: Outdoor semi-natural conditions
Twenty-four seeds from each of 21 mother trees (see Additional file 1: Table S1 in Supplementary Information) were weighed and germinated under controlled conditions (20°C, 12 h light) using a sand substrate and transplanted to outdoor semi-natural conditions. Transplanting was done periodically by blocks, in order to keep the differences due to germination included in the block effect. The experimental site was located in Madrid (40° 27' 23.82" N, 3° 45' 6.91" W, 597 m of altitude) under average environmental conditions similar to those of the seed source stand (rainfall: 452 mm, summer precipitation: 58 mm, temperature: 14.4°C, mean temperature of the coldest month: 5.3°C, mean temperature of the warmest month: 24.8°C). Actual climatic conditions during the outdoor experiment were extremely dry but similar to those in which the previous generation (i.e. the mother trees) of the population studied underwent natural selection (period 1903–1942, data not shown). As in the previous experiment, the plants were arranged in four blocks and a row-column design (within blocks), with six experimental units (one-seedling plot) for each genetic entry (24 seedlings per family). We used two contrasting planting densities; two blocks were planted at a 10x10 cm spacing (D1, 100 plants/m2), and the other two at 5x5 cm (D2, 400 plants/m2), to elicit contrasting inter-family competition effects. Seedlings were protected to avoid predation by animals (rodents, birds, or ants), and no special management was applied during the experiment duration. The experiment lasted 203 days, from April 12th to November 1st.
Phenotype measurements in seedlings
Phenotypic traits were chosen for their functional and adaptive importance in pines under limiting resource supply: growth and biomass allocation [36]-[38], specific leaf area [39] and ontogenetic change [40]. These traits present large differences among maritime pine families and populations evaluated in controlled conditions [41]-[43]. In addition, most of them seemed to have undergone adaptive genetic differentiation (as evaluated by Q
st
) in response to local climatic conditions, when tested in rangewide experiments [44].
For the two experiments, total height and shoot developmental status were evaluated in all seedlings every second week, the latter following a subjective scale based on the type of axillary shoots, as described in detail in [42]. The height of the plant at ontogenetic stage 3 (HON3) °Corresponding to the earliest observation of axillary long shoots– was used for subsequent analyses as a size-independent quantitative index of seedling development. Higher values of HON3 indicate a slower seedling ontogenetic development. Height growth of each seedling was decomposed by adjusting a Gompertz model [29] and analyzing the three main parameters of this curve: HTOT, b, and m. HTOT is the total height at the end of the experiment (asymptote), the b-parameter is related to the maximum growth rate standardized by HTOT, and m, the inflection point, indicates whether maximum growth is reached at the beginning (low values of m) or the end (high values of m) of the growing period.
At the end of each experiment, survival (SUR) was recorded, all plants were harvested, and dry mass of each fraction (needles, stem, and roots) was determined with 0.001 g precision. The dry mass of each fraction was expressed relative to that of the total plant (leaf mass fraction, LMF; stem mass fraction, SMF; root mass fraction, RMF) [36]. Total dry mass (DW) was obtained as the sum of leaf, stem, and root dry mass. Specific leaf area (SLA) was estimated from 10 needles randomly chosen from each plant. Finally, since at the seedling stage every metamer bears a primary needle, the number of stem units (NSU) was estimated by comparing total primary needle dry mass to dry mass of the 10 needles collected for SLA assessment.
Quantitative genetic models in progeny tests
For each variable and experiment, a mixed model was adjusted using Restricted Maximum Likelihood (REML):
(1)
where y
ijkl
is the phenotypic value of the variable for the lth tree from the ith family in the jth treatment located in the kth block, μ is the overall mean, ϕ
i
is the effect of the ith-family (i varying from 1 to 18–21, depending on the experiment), T
j
is the effect of the jth treatment (two watering levels in the growth chamber experiment and two density levels in the outdoor experiment, thus j varying from 1 to 2 in each experiment), χ
ij
is the family by treatment interaction, B
k
is the effect of the kth-block (1–4), c is the covariate effect, and ε
ijkl
the residual. For survival, a generalized logistic mixed model (GLMM) was applied with a binomial distribution and a logit link function. For biomass traits (LMF, SMF, RMF), the total dry mass (DW) was used as a covariate to correct for allometric changes [36]. For the rest of the quantitative traits, seed weight (WSEED) was used as a covariate to control for maternal effects. Family and treatment by family interaction were considered random factors. The models were fitted using ASREML [45].
Genetic control
Narrow-sense heritability (h2) was estimated for each trait and progeny test following standard procedures [46]. For survival (SUR), the residual variance of the GLMM was set to π2/3 [47]. The standard error of heritability was estimated using a Taylor’s series approximation [45]. Best Linear Unbiased Predictors (BLUPs) are weighted means commonly used as point estimates of random effects in mixed models. They predict the expected phenotype of a tree’s offspring or an individual by using phenotypic information collected from relatives, which is useful to evaluate genetic changes in the populations after selection. Mother-tree breeding values (based on offspring performance) for each trait were computed from the heritability of the trait and the BLUPs estimated in the experiments [45],[46]. Evolvability, as estimated by the coefficient of additive genetic variation (CVA), was also computed [10].
Phenotypic selection gradients and correlation of female reproductive success with offspring performance
Phenotypic selection gradients were computed based on traits measured directly on the mother-trees in the field, to evaluate the relationship among trait values in adults and their fitness (in terms of relative reproductive success), as obtained from molecular marker-based parentage analysis (see [26] and above). The proxy for fitness of parent trees used in this study (i.e. the actual number of successful naturally-regenerated offspring) is probably more accurate than amount of pollen donated, flower production or fruit set (e.g. [8],[13]), as the latter three do not consider post-dispersal mortality and early selection. These selection gradients were computed as the vector of partial regression coefficients of individual relative reproductive success on trait values. Given the reduced number of mother trees in the study (22), we performed a Principal Component analysis to reduce the number of dimensions in the analysis. The following log-linear model was used ([7] and references therein) to define female relative reproductive success (λ
k
):
(2)
The subscripts (1 to 3) refer to each of the three first principal components (explaining 69.64% of the variance, see Additional file 1: Table S4 in Supplementary Information). In biological terms, the β
i
are analogous to linear selection gradients, and the γ
i
are analogous to quadratic selection gradients as defined by Lande and Arnold [6]. Traits were standardized, and individual fitness was converted to relative fitness by dividing by the mean fitness of all individuals within the subset of data. The statistical significance of each specific selection gradient, β
i
or γ
i
, was estimated using likelihood-ratio tests, by subtracting the log-likelihood for the model excluding each parameter, one at a time, from the log-likelihood of the full model, which is asymptotically chi-square distributed with one degree of freedom.
Finally, to know to what extent mother reproductive success (and the associated phenotypes) is driving population phenotypic change across generations, we computed the relationship between female reproductive success and breeding values for each trait, as obtained from offspring performance under the two experimental environments, using log-linear models. Insights into this relationship are especially relevant at early establishment stages, when natural selection is stronger.
Genetic selection gradients and response to selection at early stages of establishment
Genetic selection gradients at seedling stage were computed in the two experiments using seedling height (HTOT), which is related to competitive ability and crucial for natural regeneration establishment in pines [48], and/or survival (SURV) as fitness proxies. Both fitness proxies were used in the outdoors semi-natural experiment while only seedling height was used to gauge fitness under the optimal growing conditions of the controlled-environment experiment, as no significant genetic differences among families were found for survival in this experiment (see Results). Genetic selection gradients were obtained by bivariate analysis (to reduce the bias in the estimation, see [9],[12]), as the genetic covariance between the trait and fitness divided by the genetic variance of the trait, for the traits estimated in the two experiments. In addition, we computed the genetic response to selection (R), following the Price Theorem [5], as the covariance between the trait and fitness. Standard errors were computed using a Taylor’s series approximation [45].