Functional distance is driven more strongly by environmental factors than by genetic relatedness in Juniperus thurifera L. expanding forest stands

Abstract

Key message

In our study of colonisation of abandoned farmland by Juniperus thurifera L., environmental and geographic distances were stronger drivers of functional differentiation among trees than genetic distance and tree intrinsic features. However, these effects varied both along the colonisation gradient and across predictors’ values, highlighting the complexity of the conditions under which functional diversity of forests is maximised.

Context

Forest expansion after land abandonment is changing Mediterranean landscapes, but the underlying ecological and evolutionary processes remain poorly understood.

Aims

We evaluated whether forest expansion reduces within-stand tree functional distances and whether its environmental, geographical and genetic drivers vary along the expansion gradient.

Methods

We sampled Juniperus thurifera in mature, transition and expanding stands along a colonisation gradient in Spain. We measured functional traits, genetic relatedness, surrounding vegetation cover and intrinsic features (size, phytosanitary status) on trees. We calculated within-stand pairwise functional distances and quantified the effects of genetic relatedness, geographic distance, vegetation cover and intrinsic tree features.

Results

Functional and genetic distances showed no consistent patterns along the expansion gradient. Functional distance was unrelated with genetic distance and strongly driven by geographic distance, likely representing unmeasured environmental factors. Vegetation cover and intrinsic tree features affected functional traits in transition and expanding stands. The effects on functional distance changed across predictors’ values.

Conclusion

The role of environmental and geographic factors on tree functional divergence varied along the farmland colonisation gradient, while we did not detect effects of genetic distance. Accounting for this context-dependence and for non-linear effects on functional traits can improve our understanding and management of forest expansions.

1 Introduction

European countries are experiencing important changes in land uses associated with rural exodus (Rudel et al. 2005; Song et al. 2018). The abandonment of cultivated land is leading to a widespread expansion of second-growth forests, a process with particular relevance for landscape and ecosystem functioning in Mediterranean regions (Valladares et al. 2014; Doblas-Miranda et al. 2017; Palmero-Iniesta et al. 2020). Forest expansion is of high relevance due to the diversity of ecosystem services provided by forests to human populations, including timber and non-timber products as well as the regulation of climate, pests, floods and nutrient cycling (Gamfeldt et al. 2013; Navarro and Pereira 2015; Vilà-Cabrera et al. 2017). However, we are far from understanding the ecological and evolutionary processes that underlie the expansion of second-growth forests and how they influence the functional and genetic properties of these systems.

Analysing functional diversity can be a tool for understanding and predicting the consequences of human activities on ecosystem functioning, such as changes associated to land abandonment (Lavorel et al. 2011; Wellstein et al. 2011). Functional traits (i.e. morphological or physiological features that influence individual fitness; Reich et al. 2003; Violle et al. 2007) have been used for assessing community assembly (McGill et al. 2006), and the success of forest restoration (Martínez-Garza et al. 2013). At the intraspecific level, variability in functional traits defines a species’ capacity to optimise fitness along environmental gradients and, ultimately, the breadth of its ecological niche (Ackerly and Cornwell 2007; Albert et al. 2010). Intraspecific trait variation informs thus about the ecological success of populations and species and their ability to persist upon environmental changes (Wang et al. 2010; Forsman and Wennersten 2016). In second-growth forests, studying the effects of expansion into abandoned land on functional traits can also be crucial for assessing species vulnerability to extreme climatic events (Alfaro-Sánchez et al. 2019). Overall, disentangling the patterns and drivers of functional diversity in the tree species leading forest expansion is important for a better management of the landscape and ecosystem services.

Functional diversity itself but also its intrinsic and extrinsic drivers are likely to vary along the forest expansion process. The colonisation of new areas typically occurs through relatively few pioneer individuals that arrive by long-distance dispersal and fill the neighbourhood with their progeny (Petit et al. 2004; Hampe 2011). As a consequence, genetic diversity may be reduced at the colonisation front owing to founder effects (Excoffier et al. 2009). A high genetic relatedness of the individuals near the colonisation front could in turn lead to high functional similarity, resulting in low stand-level diversity in functional traits. Diverse ecological, physiological and morphological traits have shown phylogenetic conservatism at the intraspecific level (Wiens et al. 2010; Pearman et al. 2010). In forest trees, numerous common-garden experiments and provenance trials have also documented that functional traits are to some extent genetically determined, although this effect is often moderate owing to the polygenic inheritance and phenotypic plasticity of many traits (Alberto et al. 2013; Benito Garzón et al. 2019).

In more mature forests located behind the colonisation front, functional traits may show lower genetic effects and stronger signals of acclimation to surrounding biotic and abiotic conditions, e.g. through phenotypic plasticity (Niinemets 2010). Surrounding neighbours of the same or different species may be particularly influential for tree establishment (Magee and Antos 1992; Rice et al. 2012) and functional traits (Jose et al. 2004; Coop and Givnish 2008). Functional traits might also be conditioned through physiological or allometric constraints by intrinsic characteristics of the trees, such as their size (Meinzer et al. 2011; Hietz et al. 2017) and health status (Gilbert 2002; Tulik 2014). Yet, the importance of these intrinsic features on tree functionality may decline once acclimation to surrounding environment has taken place in mature forest stages. Neighbour competition seems indeed more determinant for tree performance than individual size in Mediterranean forests (Gómez-Aparicio et al. 2011).

Here we assessed the patterns and drivers of functional traits in the Mediterranean tree Juniperus thurifera, along a gradient of colonisation. Juniperus thurifera forests have been vigorously expanding in central Spain during the last decades as a consequence of widespread land abandonment (Blanco et al. 2005; Olano et al. 2008). We established multiple study plots in mature, transition and expanding forest stands and measured leaf and wood functional traits, molecular genetic markers, environmental conditions (vegetation, conspecific and rock cover) and intrinsic tree features (size and phytosanitary status). The first objective of our study was to quantify functional and genetic distances between the individual trees within each plot, hypothesising that these distances would be lowest in the expanding forest stage. The second objective was to analyse the potential drivers of functional distance between individuals and whether these drivers changed along the forest expansion gradient. We predicted that functional distance would be mainly related to (1) genetic distance and distance in intrinsic tree features in the expanding front and (2) environmental distance in the mature forest. We also analysed tree age across forest stands for a better interpretation of results.

2 Material and methods

2.1 Study species and area

Juniperus thurifera L. (Cupressaceae) is a wind-pollinated dioecious tree with individuals growing 5–10(− 20)-m tall (Montesinos 2007). The species forms low-density forests on poor, shallow, rocky soils that can establish at high altitudes and cover broad climate ranges (Gauquelin et al. 1999). Juniperus thurifera tolerates broad temperature ranges and is adapted to water stress typical of Mediterranean-climate regions, where precipitation is moderate or low and it is concentrated in spring (Montesinos et al. 2009). Its extant distribution area is restricted to the Western Mediterranean Basin with populations in Spain, Morocco, France, Algeria and the Italian Alps. Approximately 90% of the species distribution is located in Spain, where it occupies 600,000 ha with 117,000 ha being monospecific forests (Alonso Ponce et al. 2010). The present study was carried out in the Alto Tajo Natural Park and surrounding areas (Guadalajara province, central Spain). The climate is continental Mediterranean with hot dry summers and cold snowy winters (mean annual rainfall ± SE, 477.1 ± 15.6 mm; mean annual temperature, 10.4 ± 0.2 °C; Molina de Aragón 40°50'40" N, 1°53'07" W, 1063 m.a.s.l., 1951–2017 period; Spanish National Meteorology Agency). The effective population size of J. thurifera in the study area has not been quantified, but the adult population is extensive and relatively continuous, and its effective size can therefore be considered to be large (i.e. thousands of individuals). In our study area, J. thurifera seeds are dispersed by a rich community of vertebrates, including thrushes (dispersal distances up to a few hundred meters, conditional on the presence of tree perches) and carnivores (dispersal distances up to 1–2 km, irrespective of tree presence; Escribano-Avila et al. 2014).

2.2 Field sampling

We established 17 plots representative of a gradient of forest expansion after land abandonment, located in three study sites (Maranchón, Huertahernando and Ribarredonda; Table 1; Fig. 1). Five plots were located in mature forest stands (a core of well-conserved forest), six plots near the front of expansion (on abandoned agricultural lands) and six plots in a transition zone between the two former stages. Mean plot size varied between mature (0.51 ha), transition (0.81 ha) and expanding (1.35 ha) stands to account for differences in density of J. thurifera adults. Mean distance among plots was 1.5 km in Maranchon, 0.4 km in Huertahernando and 0.3 km in Ribarredonda (due to differences in site accessibility). We conducted field sampling between September and October 2017. In total, we randomly selected 394 adult individuals across study plots (23 individuals per plot on average; Table 1). We considered adult those trees higher than 1.40 m and with a trunk diameter at breast height (DBH; 1.3 m) equal or greater than 3 cm (the minimum size for which we found reproductive individuals). For each tree, we took the geographic position (latitude and longitude) and measured four biometric variables: tree height (measured with Haglöf Vertex IV hypsometer), crown projection (calculated as the mean of the projection of two perpendicular axes that passed through the axis of the trunk; measured with a Haglöf DME distance measurer), mean DBH across all stems and quadratic mean diameter (Curtis and Marshall 2000). We cored each individual twice at 50 cm from the ground using a Pressler increment borer (0.5 cm; Haglöf, Långsele, Sweden) to obtain wood cores with intact pith that were brought to the laboratory for the determination of tree age and wood density. Moreover, we visually estimated the phytosanitary status of the tree using a semiquantitative scale with three categories: (1) individual with several dry branches; (2) individual with one dry branch and (3) individual with no dry branches or apparent damage. For each individual, we delimited the area of influence as a circle around it with a radius equal to crown projection (see above). Within this area, we estimated the percentage of cover occupied by individuals of other woody species (‘vegetation cover’), individuals of J. thurifera (‘conspecific cover’) and rocks. This was made for the north-facing and south-facing halves of the area of influence, and the final estimation for each type of cover was the average of north and south values. Vegetation and rock cover were strongly correlated in a linear model (slope = − 0.88, p < 0.001, R² = 0.71), so we discarded percentage of rock cover from the posterior analyses.

Table 1 Study plots of Juniperus thurifera in three stages along the forest expansion gradient (mature, transition and expanding) in Guadalajara, Spain

Fig. 1

Location and pictures of Juniperus thurifera study area. a We selected three sites in central Spain (Guadalajara, in purple). Study plots are located in mature, transition and expanding forest stages in b Maranchón, c Huertahernando and d Ribarredonda (see also Table 1). Study plots differ in J. thurifera density between e mature, f transition and g expanding stages (pictures correspond to Huertahernando)

For the analyses of leaf functional traits, we collected 20 foliar units formed during the previous year from each individual (Villellas et al. 2020a). All foliar units came from branches oriented towards the south and located in the upper part of the canopy to minimise unaccounted micro-environmental variability. For the genetic analyses, we collected apical leaves from all the 394 individuals sampled for functional traits, plus an additional set of 269 georeferenced adult individuals from the same study plots to increase sample size (in total, we collected samples for genetic analyses from an average of 39 individuals per plot; Table 1).

2.3 Sample processing

To estimate tree age, one core per tree was air dried, glued onto wooden mounts and polished using sandpaper of progressively finer grain until tree rings were visible. Tree rings were counted with a stereomicroscope; cores were scanned at 1600 d.p.i. and ring widths measured to an accuracy of 0.001 mm using the software CooRecorder v9.3 (Cybis Elektronik & Data AB, Saltsjöbaden, Sweden). Cross-dating of individual series was checked using the software CooRecorder and COFECHA (Holmes 1983). The other core from each individual was used to calculate wood density as dry weight per unit of saturated volume. We measured dry weight after oven-drying cores at 60 °C for 72 hours. Then, we rehydrated cores by introducing them in distilled water for at least 24 hours. Finally, we calculated saturated volume by introducing cores in a container filled with water and by measuring the weight of displaced water with a precision balance.

We measured a set of leaf traits with key ecophysiological function: leaf dry matter content (LDMC), carbon isotope signature (δ¹³C), carbon (C) and nitrogen (N) concentration and carbon:nitrogen ratio (C:N). LDMC and C concentration inform about leaf density and structural support, δ¹³C is related to water use efficiency, N concentration informs about photosynthetic capacity and C:N ratio is related to nitrogen use efficiency (Dawson et al. 2002; Cornelissen et al. 2003; Dijkstra et al. 2016; Carvalho et al. 2020). To measure LDMC, we pooled ten foliar units collected from each individual. Leaves were rehydrated during 12 h to obtain their saturated weight and oven-dried them at 60 °C for 72 h to obtain the dry weight. To calculate LDMC, we divided saturated weight by dry weight. The other ten foliar units were finely ground to be used for isotope analyses (MM300, Retsch Haan, Germany). We weighed a subsample of 3–4 mg and placed it into a tin capsule for the analysis of δ¹³C (‰) and leaf C and N concentration (μg mg⁻¹sample). Analyses were carried out at the UC Davis Stable Isotope Facility using a PDZ Europa ANCA-GSL elemental analyser interfaced to a PDZ Europa 20-20 isotope ratio mass spectrometer (Sercon Ltd., Cheshire, UK). Carbon isotope signature (δ¹³C, ‰) values were defined as (R_sample/R_standard) × 1000, where R_sample was ¹³C/¹²C ratio for leaf samples and R_standard was international standard ¹³C/¹²C ratio (Vienna Pee Dee Belemnite, VPDB). Finally, we divided leaf C concentration by N concentration to calculate C:N ratio.

2.4 Genetic analyses

We used the ddRAD-Seq genomic data described by García et al. (2018) to design six multiplexes using MassARRAY Assay Designer v 4.0.20.2 (Agena Biosciences, USA), for a total of 239 SNPs. This screening step was performed on a random subset of 95 samples using the iPLEX Gold chemistry on the MassARRAY System of the Genome Transciptome Facility of Bordeaux. Data analysis was achieved using MassARRAY Typer Analyser 4.0.20.65. We filtered out all loci that were either monomorphic, showed a weak or ambiguous signal (i.e. with unclear cluster delimitation or displaying more than three clusters of genotypes) or contained too many missing data. Two final multiplexes of 38 SNP each were then designed using the best markers from the screening step and were used for the genotyping of all individuals (Villellas et al. 2020b). The corresponding data analysis was also performed with MassARRAY Typer Analyser. After establishing the genotypes, we removed all loci and all individuals with > 20% of missing data. The resulting final dataset contained 649 individuals genotyped at 55 loci and contained 3.15% of missing data.

Juniperus thurifera is a tetraploid species, and it was impossible to reliably identify the different allele combinations of each locus and their respective dosage. Therefore, we coded all genotypes instead as incomplete polyploid individuals for the following analyses (i.e. we retained three classes for each locus: two corresponding to the homozygotes of each allele and the third pooling all heterozygotes). Nason’s kinship coefficient (Loiselle et al. 1995), a commonly used estimator of relatedness between individuals that assumes codominance and does not weigh allele frequencies, was calculated with the programme SpaGeDi version 1.2 (Hardy and Vekemans 2002). Random permutation of genes among all individuals and random permutations of individual locations among all individuals were performed to test for the significance of kinship coefficients. For this analysis, we only considered one distance class (all tree pairs located within 10 m) to account for differences in density among plots. In turn, we pooled the plots for a given stage and study site within the same kinship analysis.

2.5 Statistical analyses

2.5.1 Cover and age analyses

We evaluated whether plots differed between stages of the forest expansion gradient in terms of conspecific cover and tree age. To this aim, we used linear mixed models (LMM) with stage and site as fixed factors and study plot as a random effect (see Table 2 in the Appendix for equation). We applied a square-root transformation to conspecific cover to normalise residuals. For both analyses, we finally performed Tukey tests among stages.

2.5.2 Distance analyses

Our principal statistical analyses relied on the distances—functional, genetic, geographic, environmental and intrinsic—between all pairs of individual trees growing within a given study plot. We first quantified Euclidean distances in functional traits, both for each trait separately (C:N ratio, δ¹³C, N, C, LDMC and wood density) and for all traits together. Second, we quantified their genetic distances multiplying pairwise kinship coefficients by − 1. Third, we computed trees’ geographic distances. Finally, we calculated Euclidean distances for vegetation cover, conspecific cover, phytosanitary status and tree size. The variable tree size was derived conducting a principal components analysis on our four biometric variables (mean DBH, quadratic mean diameter, crown projection and height); we used the first principal component that explained 82.4% of the overall variance as a compound measure of tree size. To evaluate how distances in functional traits and their potential predictors (genetics, geographic position, environment and intrinsic tree features) varied along the forest expansion gradient, we performed LMMs including forest stage and site as fixed effects and study plot as a random effect (see Table 2 in the Appendix for equation). We then tested for differences among stages with Tukey contrasts and calculated mean predicted values per forest stage. These analyses were only performed for the individuals with functional trait data.

To analyse the drivers of functional distance between individuals, we used Generalised Dissimilarity Models (GDM). GDM is an extension of matrix regression based on linear combinations of I-spline basis functions, which accommodate non-linear relationships between predictor and response distance matrices (Ferrier et al. 2007). The resulting monotonic functions (one for each predictor) have two key properties: (1) their maximum height represents the change in the response variable explained by each predictor; and (2) their slope indicates the rate of change in the response variable along the predictor’s range of values (Ferrier et al. 2007). GDMs also quantify the effects from geographical distances between individuals not accounted for by other predictors included in the models. For each plot, we performed a GDM to test whether distances between individuals in functional traits were explained by genetic, environmental and/or geographic distances. We also included distances in tree size and phytosanitary status as explanatory factors in GDM. All distances were scaled to a range from 0 to 1 before the analyses. In each model, the significance of each predictor was evaluated with 250 permutations of its corresponding distance matrix. Permutations allowed us to compare the observed difference in explained deviance between models with and without the predictor to a null distribution of differences in deviance (obtained with the permuted distance matrices). The significance of predictors was evaluated in models that retained only the most important factors. For that aim, we used a backward-elimination method where each GDM was repeated without the least influential factor until all predictors present in the model had a p value ≤ 0.2 (or until the model had a minimum of two predictors even if their p values were higher than 0.2). We verified that changing this threshold p value across a range of 0.1–0.3 did not modify results. In 14 out of 119 GDM performed in total, the model did not converge for one of the explanatory factors due to a weak relationship between the response and predictor. For each functional trait and forest stage, we evaluated the number of plots in which a predictor showed a significant (p < 0.05) or marginally significant (0.05 < p < 0.1) effect. All analyses were performed with the packages nlme, multcomp, gmt, lsmeans and gdm in R software version 3.4.0 (R Development Core Team 2017).

3 Results

3.1 Cover and age analyses

Mean percentage cover by J. thurifera adult trees was 20% in the mature stage, 10% in the transition stage and 7% in the expanding stage. Cover by J. thurifera differed between mature and transition stages (z = − 3.58, p < 0.001), between mature and early stages (z = − 5.89, p < 0.001) and between transition and early stages (z = − 2.34, p = 0.049). Mean tree age was 44 years in the mature stage (range, 13–116), 30 years in the transition stage (range, 12–78) and 28 years in the expanding stage (range, 11–88; Fig. 5 in the Appendix). Tree age differed between mature and transition stages (z = − 2.72, p = 0.018) and between mature and early stages (z = − 2.86, p = 0.012), but not between transition and early stages (z = − 0.10, p = 0.994).

3.2 Distance analyses

Trees tended to differ more in the mature stage than in other stages for most functional traits, but none of these differences was statistically significant (Fig. 2a–g). Individuals were more distant geographically, and more different in phytosanitary status, in the expanding stage than in the transition or mature stages; other predictors (genetic, cover and size distances) showed no significant differences among stages (Fig. 2h–m). For the subset of pairs of individuals located closer than 10 m, permutation tests showed diverging results among sites: the highest (and significant) mean kinship coefficients corresponded to the expanding stage in Huertahernando and to the mature stage in Ribarredonda, whereas no significant patterns were found in Maranchón (Table 3 in the Appendix).

Fig. 2

Mean distances between individuals of J. thurifera in functional traits (a–g) and in different types of predictors (h–m). Results are mean predicted values (with 95% confidence intervals) for the three stages of forest expansion: mature (Mat), transition (Trn) and expanding (Exp). Functional traits shown are carbon isotope signature (δ¹³C), the ratio between leaf carbon and nitrogen concentrations (C:N), leaf carbon concentration (C_con), leaf nitrogen concentration (N_con), leaf dry matter content (LDMC), wood density (W.DEN) and a combination of the individual traits studied (COMB). Predictors are genetic distance (GEN), geographic distance (GEO), distance in vegetation cover (VEG.COV), distance in conspecific cover (CON.COV), distance in tree size (SIZE) and distance in phytosanitary status (PHYT.ST). Values are Euclidean distances except for genetic distance, which corresponds to kinship multiplied by − 1, and geographic distance, which is shown in meters. Lowercase letters in I) and M) denote significant differences in distance among forest stages (stages not sharing a letter are significantly different). Note differences among panels in vertical axis scale

Functional distance among trees was affected by different predictors depending on the stage of the expansion gradient (Fig. 3). In most traits, the GDMs retaining the highest number of significant or marginally significant predictors corresponded to plots in the transition stage (Fig. 3; Table 4 in the Appendix). Genetic distances had virtually no effects on distances in functional traits. Geographic distance showed significant or marginally significant effects in two or more plots on all functional traits in the transition stage and on several traits (C:N ratio, δ¹³C, N concentration and wood density) in the mature and expanding stages. Vegetation cover was a frequent predictor for some traits (C:N, δ¹³C, N concentration and the combined trait) but only in the transition stage. Conspecific cover, tree size and phytosanitary status affected several traits, but the effects in more than one plot were restricted to the younger forest stages: in the expanding front, conspecific cover affected δ¹³C and C concentration, tree size affected δ¹³C and phytosanitary status affected C concentration; in the transition stage, conspecific cover affected C:N and LDMC, tree size affected C concentration and phytosanitary status affected C:N, wood density and the combined trait. The GDMs performed with the combination of all functional traits did not show qualitatively different results than GDMs performed with single traits (Fig. 3g): all the predictors appearing as frequent for the combined trait were also frequent for other traits, in particular for C:N (the trait with the highest number of significant or marginally significant effects of predictors on functional distance).

Fig. 3

Significance of predictors of functional divergence in Juniperus thurifera, represented with p values from Generalised Dissimilarity Models (GDM). Functional traits (rows) and predictors (columns) are as in Fig. 2. For each functional trait, a GDM including all predictors was performed in each study plot of mature (Mat), transition (Trn) and expanding (Exp) forest stages. The number of plots per forest stage is 5 (Mat), 6 (Trn) and 6 (Exp), but p values of predictors are only shown for plots with significant (dark blue dots; p < 0.05) or marginally significant (light blue dots; 0.05 < p < 0.1) effects. For each forest stage, predictors are shadowed in grey when they show significant or marginally significant p values in two or more plots. Vertical axes are in square-root scale to improve display

The I-spline functions of the GDMs showed that predictors differed in their importance for functional traits and in the shape of their curves (Fig. 4; Figs. 6, 7 and 8 in the Appendix). Despite the high number of GDMs that revealed significant effects of geographic distance on functional distance, the effects of this predictor were generally less important than the effects of other significant predictors. In contrast, distance in vegetation cover was not a frequent predictor of functional distance but, when present, showed the highest importance. The importance of other predictors (conspecific cover, tree size and phytosanitary status) varied among GDMs. The shape of I-spline functions of predictors showed that the largest effects on functional distance (i.e. the strongest slopes) tended to occur at low values of tree size, at high values of geographic distance and all along the range of values of conspecific cover, although there were discrepancies to these generalities. The functions of vegetation cover and phytosanitary status showed no clear patterns in slope along the range of predictor values.

Fig. 4

Results from Generalised Dissimilarity Models (GDM) for two of the studied functional traits: carbon isotope signature (δ¹³C; a–c) and the ratio between leaf carbon and nitrogen concentrations (C:N; d–f), and for each functional trait, rows correspond to mature (Mat), transition (Trn) and expanding (Exp) forest stages, and panels correspond to GDMs performed in each plot (with plot name indicated at the top of the panel). Curves represent the I-spline functions of predictors with significant (solid line; p < 0.05) or marginally significant (dotted line; 0.05 < p < 0.1) effects on divergence in functional traits between individuals. Predictors are as in Fig. 2 (genetic distance is not included due to a lack of significant effects in the selected functional traits). The values of predictors in horizontal axes were scaled to a 0–1 range to show effects from different predictors in the same panel. The maximum height of the I-spline functions can be interpreted as the importance of the predictor, and the slope provides an idea of the range of values of the predictor in which the effect on functional divergence is higher

4 Discussion

Disentangling the patterns and drivers of functional diversity in forest expansion is of great importance for understanding and predicting changes in Mediterranean landscapes (Valladares et al. 2014; Doblas-Miranda et al. 2017). In our study of expanding J. thurifera forests, we found no clear differences in functional or genetic distances between individuals among the different stages of the expansion gradient. Functional divergence was not related to genetic distances and was mainly driven by geographic distance between individuals instead. Differences in surrounding vegetation, tree size and phytosanitary status also promoted divergence in some functional traits but only in the transition and expanding forest stages.

In contrast to our expectations, individuals in the expanding and transition stages were on average younger that in the mature forest, but they did not show consistently lower genetic distances (higher kinship). In fact, only one of the three sites showed higher kinship in the expanding front than in the mature stand. Our results suggest that the colonisation process occurs with enough immigrant gene flow to preclude any bottlenecks, although additional analyses would be needed for confirmation. Studies in leading-edge populations of wind-pollinated trees suggest that genetic diversity can be rapidly restored through the interacting effects of long-distance pollen flow and purging of inbred individuals during recruitment (Hampe et al. 2013). Another mechanism that could have minimised the effects of founder events in our study sites is effective long-distance seed dispersal from mature forest stands in the surrounding landscape matrix (Petit et al. 2004; Jordano et al. 2007; Pesendorfer et al. 2016). Previous work in J. thurifera expanding stands has indeed shown a rich community of seed dispersers including species known to perform regular long-distance dispersal (Escribano-Avila et al. 2012, 2014,), as well as effective plant recruitment away from reproductive adults (Gimeno et al. 2012).

The apparent absence of founder effects may explain the lack of relationships between genetic and functional distances. Any eventual effects of genetic relatedness seem in fact overshadowed by trait acclimation to local conditions. Differences between genetic and functional patterns may also arise because genetic variation was characterised with neutral markers; functional traits, unlike neutral genetic markers, are in fact driven not only by random processes but also by selection and multiple gene interactions (van Tienderen et al. 2002; Dlugosch and Parker 2008). However, it should be noted though that our failure to detect links between genome-wide relatedness and phenotypic trait variation does not necessarily imply that the latter is not under differential selection, because our limited marker coverage does not ensure linkage to causal polymorphisms underlying trait variation. We are not aware of similar studies in trees, but some recent work in herbaceous and shrub species has shown stronger effects of epigenetic rather than genetic variation on functional traits (Herrera et al. 2017; Lele et al. 2018). The role of epigenetic variation during tree expansion deserves thus further attention, particularly in relation with phenotypic plasticity (Zenni et al. 2017).

The present study confirmed our previous expectation that patterns of functional divergence in the expanding front partly responded to intrinsic characteristics of the trees. The influence of tree size and health status has been frequently reported in previous studies of tree performance (Gilbert 2002; Meinzer et al. 2011; Tulik 2014; Hietz et al. 2017) and should be considered when evaluating functional diversity patterns. On the other hand, neighbour presence significantly influenced functional traits in the transition and expanding forest stages, in contrast with our preliminary hypothesis. Previous research in J. thurifera detected both facilitative and competitive effects of conspecific cover on tree growth, particularly in saplings and juveniles (Gimeno et al. 2012, 2015). In our study, we also found effects of conspecific trees on adult individuals, particularly on functional traits like C:N, δ¹³C and C concentration. Moreover, we detected an effect of other woody species on some traits. Notably, the effect found on leaf N concentration and C:N ratio could be due to soil fertilization by co-occurring nitrogen-fixing shrubs such as Genista species (Vidal et al. 2019). Thus, it seems that biotic interactions with both conspecifics and individuals of other species are driving patterns in J. thurifera functional traits. However, we cannot discard that some of the observed patterns resulted from independent responses of coexisting individuals and species to the same abiotic conditions.

While surrounding vegetation and intrinsic characteristic of trees were important in the younger forest stages for some functional traits, the main driver of pairwise functional divergence in J. thurifera stands was geographical distance. Our results confirm the importance of considering spatial autocorrelation in the analyses of functional diversity (Legendre 1993; Biswas et al. 2016). Interestingly, several traits showed significant spatial autocorrelation in the mature forest stage, where no other factor exerted significant influence on functional divergence. The effect of geographical proximity probably reflects the consequences of additional biotic or abiotic factors unmeasured in our study. For example, soil or topographic conditions, which frequently show local-scale variation, have been considered determinant for functional traits in trees (Fayolle et al. 2012; Schneider et al. 2017). However, whether spatial autocorrelation in the older and younger stands of J. thurifera respond to the same underlying processes remains an open question. Further studies would be needed to ascertain the environmental drivers of functional traits not yet identified and also to clarify the role of phenotypic plasticity vs. genetic differentiation in the observed trait-environment relationships.

Two additional results from our analyses deserve special attention. First, most effects detected in our study were stronger in the transition than in the expanding and mature forest stages. In the expanding front, which consists of young forest stands in development, the influence of random factors on functional traits could be more important than in other stages of the gradient. In contrast, traits in the mature forest could be driven by acclimation to a few but highly influent local environmental factors, which would dominate over other potential sources of functional variation. These results highlight the importance of studying ecological and evolutionary processes along different steps of the expansion gradient rather than using binary comparisons (e.g. mature vs. expanding). Spatial and temporal shifts in the main drivers of forest expansion have also been reported in previous studies (Hansen et al. 2016; Guariguata and Ostertag 2001). In second place, we found that the effects on functional divergence varied along predictors’ range of values, with the largest effects taking place in general at low values of tree size and high values of geographic distance. The importance of non-linear relationships for the ecological and physiological processes underlying forest growth and expansion has been previously highlighted (Burkett et al. 2005). Overall, these comprehensive analyses allow us to identify not only the major drivers of variation in functional traits but also the stages along colonisation gradients and the range of predictors’ values maximising functionally diverse forest stands. These are valuable considerations for conservation and restoration practices, in view of the importance of intraspecific trait variation for species niche breadth and their ability to persist upon environmental changes (Ackerly and Cornwell 2007; Forsman and Wennersten 2016).

5 Conclusion

Our study of colonisation of abandoned agricultural lands in the Mediterranean tree J. thurifera detected no consistent genetic differences along the forest expansion gradient. Wind pollination and a rich community of fruit dispersers have likely contributed to a continuous arrival of new genes to the expansion front from nearby populations. The markers we used did not allow us to detect relationships with functional divergence among trees. Tree functional distance was instead mainly determined by geographic location, suggesting the effects of unmeasured local environmental factors. The drivers of functional divergence changed along the expansion gradient, with neighbour presence and intrinsic tree features having a stronger influence in younger than in more mature forest stands. Finally, we showed that accounting for non-linear patterns in functional traits can be very useful for a better understanding and management of functional diversity in a context of forest expansion after land-use changes.

Appendix

Table 2 Equation of the models performed to test for differences among mature, transition and expanding forest stages for a set of response variables: conspecific cover, tree age, distances in functional traits (carbon isotope signature, the ratio between leaf carbon and nitrogen concentrations, leaf carbon concentration, leaf nitrogen concentration, leaf dry matter content, wood density and a combination of the individual traits) and distances in predictors (genetic distance, geographic distance, distance in vegetation cover, distance in conspecific cover, distance in tree size and distance in phytosanitary status)

Table 3 Mean kinship coefficients (and lower and upper limits of 95% confidence intervals) for pooled plots within mature (Mat), transition (Trn) and expanding (Exp) stages of each study site of Juniperus thurifera

Table 4 Deviance explained by Generalised Dissimilarity models (GDM) for each functional trait of Juniperus thurifera in plots at three forest stages

Fig. 5

Frequency of ages of individual trees of Juniperus thurifera in plots of mature (light grey), transition (dark grey) and expanding (black) forest stages

Fig. 6

Results from Generalised Dissimilarity Models (GDM) for two of the studied functional traits leaf carbon concentration (C_con; a–c), leaf nitrogen concentration (N_con; d–f). See Fig. 4 in the main text for further details. Note that the pink dotted line corresponds to the only (marginally significant) effect found for genetic distance across all GDMs performed in the study

Fig. 7

Results from Generalised Dissimilarity Models (GDM) for two of the studied functional traits: leaf dry matter content (LDMC; a–c) and wood density (W.DEN; d–f). See Fig. 4 in the main text for further details

Fig. 8

Results from Generalised Dissimilarity Models (GDM) for the combination of all studied functional traits (COMB). See Fig. 4 in the main text for further details