Large-Scale Regeneration Patterns of Pinus nigra Subsp . salzmannii : Poor Evidence of Increasing Facilitation Across a Drought Gradient

Tree recruitment is a key process underlying stand dynamics and sustainability in managed forests. Woody plant cover is known to affect the regeneration success of Pinus nigra, suggesting the existence of facilitative plant-plant interactions. The regeneration patterns of this Mediterranean pine were analyzed across its distribution area, using data from 3226 plots of the Spanish National Forest Inventory. We aimed to test the hypothesis that seedlings establishment occurs under higher values of either canopy or shrub cover in the driest populations, as predicted by the stress-gradient hypothesis. Data were analyzed by means of Generalized Linear Models and multivariate methods. Results revealed that regeneration failure occurs on a regional scale, and that regeneration is facilitated by tree canopy cover of 55%–80%. A non-linear pattern of interaction along an aridity gradient was identified, with competition at the wettest site, high facilitation at the mid-dry sites, and low facilitation at the driest site. Evidence suggests that some shrub species may facilitate recruitment in the harsher areas. Collectively, our results reduce the possibilities of adapting forest management to drying climates by the application of alternative silvicultural prescriptions involving canopy cover.


Introduction
Water shortage has been widely reported as the main cause impeding the establishment of juvenile plants in Mediterranean ecosystems [1][2][3].Seedlings are highly vulnerable to drought stress, and Mediterranean-type ecosystems are characterized by a sharp seasonality of rainfall and temperature, with hot-dry summers [4].Consequently, seedling mortality is found to be exacerbated both during the summer season and in open areas, where high irradiance increases evaporative demands and further reduces soil-water availability [5][6][7][8].Studies focused on the recruitment dynamics of Mediterranean tree species have found increasing seedling mortality during the summer, compared with other periods of the year, and in open areas, compared with sites located under the canopy of woody plants [9,10].Thus, Mediterranean tree species generally benefit from being planted in the proximity of shrubs during reforestation work, as most shrub species may buffer against high temperatures and radiation [8].Similarly, the risks of excessive evapotranspiration, overheating, and photo-inhibition are ameliorated under the cover of tree canopies [2,10,11].These positive plant-plant interactions between tree seedlings and woody species are examples of facilitation [12].
According to the stress-gradient hypothesis [13], facilitation should increase across gradients of abiotic stress, i.e., the shade of neighboring plants could be expected to increasingly ameliorate drought stress as conditions become drier.This prediction has been validated for a broad range of environments and species, including species of the Pinus genus [14][15][16].However, the stress-gradient hypothesis has been challenged by empirical evidence coming from arid and semi-arid areas.In these dry environments, drought stress might increase in the vicinity of other plants, despite the benefits of shading, as neighboring plants reduce the overall water content of the soil by direct water uptake and rainfall interception [17,18].Furthermore, the net facilitative outcome of shrub-seedling interactions along a climate gradient will depend on the interacting shrub species, as the nurse effect provided by a shrub is modified by its morphological traits, such as size, branch disposition, and spininess [19], and its phylogeny [14,20].
Contrasting with the stress-gradient hypothesis, the trade-off hypothesis predicts that shade will increase the negative effect of drought, as water balance and carbon uptake cannot be maximized simultaneously [21].Plants living in poor light environments would need to allocate more resources to aerial parts, and less to roots, thereby increasing their potential for light interception, while diminishing their ability to capture water and survive drought.The trade-off hypothesis has been supported by different ecophysiological studies involving Mediterranean species [22,23].As a way to reconcile the apparent opposing views of the stress-gradient and trade-off hypotheses, it has recently been proposed that drought effects are ameliorated at intermediate irradiance, while they are more severe at higher or lower light levels due to a switch from facilitation to competition at both ends of an aridity gradient.When this situation applies, the expected plant response to drought along the irradiance gradient should be non-linear (a humped-back shape model) [24].However, it is also possible that biotic interactions become unimportant relative to the effect of the environment when environmental conditions become very severe, generating a collapse of facilitation [25].Whatever the hypothesis considered, an understanding of facilitative plant-plant interactions is important for forest managers, as tree harvest and/or shrub clearance modify the amount of light reaching the forest floor and, thus, forest regeneration dynamics.This understanding might be of greater importance in a climate change scenario, which implies rising temperatures and altered rainfall patterns in several Mediterranean-type forest ecosystems [26].In this respect, there is increasing concern about the long-term persistence of some managed Mediterranean mountain pinewoods, which exhibit limited regeneration dynamics [2,27].
Pinus nigra Arn.subsp.salzmannii Dunal (Franco) (P.nigra hereafter) is a Mediterranean mountain pine, commonly harvested [28], and its forests are included in the European Union list of natural habitats requiring specific conservation measures.Several small-scale studies have shown that P. nigra regeneration frequently fails under natural conditions, particularly in open microhabitats, but it is unknown whether such lack of regeneration is a common pattern at the regional scale [6,9,27].The species is widely distributed along the eastern mountain ranges of the Iberian Peninsula.Here, summer precipitation decreases and temperature increases from north to south, resulting in a latitudinal gradient of aridity [29].We hypothesize that the establishment of P. nigra seedlings should be more closely related to either canopy or shrub cover in the southern populations than in the northern ones, in accordance with predictions from the stress-gradient hypothesis.
The specific questions addressed in this study were: considering results from small-scale studies, ( 1) at what scale does regeneration failure of P. nigra occur, and ( 2) what are the biotic and abiotic factors determining P. nigra regeneration at the regional scale?(3) Is P. nigra regeneration facilitated by tree canopies and/or shrubs, and, if the latter applies, what are the patterns of facilitation across the species distribution area?

Species, Data Set, and Study Variables
P. nigra thrives along the eastern calcareous mountains of Spain (Figure 1), under a Mediterranean mountain type climate characterized by cold-wet winters and hot-dry summers [26].Seed predation, herbivory by ungulates, the presence of vegetation litter and the irregular production of seed crops are factors which negatively affect P. nigra regeneration in natural stands, although summer drought is the major "bottleneck" of P. nigra demography [27,28,30,31].
From the 3SNFI provincial databases, we selected inventory plots that (i) were situated within natural stands; (ii) contained at least one P. nigra individual with a d.b.h.≥20 cm (in order to assure a seed source); (iii) and were located within the provinces of Lleida, Teruel, Cuenca, and Jaé n-Granada.This selection of provinces allowed us to cover a latitudinal gradient of 650 km, including data from the mountainous areas with the most extensive forests of P. nigra (study regions hereafter): Prepirineo Catalán (CAT), Alto Maestrazgo-Montes Universales (MAE), Sistema Ibé rico Meridional (IBE), and Cordilleras Bé ticas (BET) (Figure 1).A total of 3226 plots were selected by this filtering process.The region, the abundance of P. nigra regeneration (seedlings and saplings), and the corresponding values of four stand-structural variables and six soil-physiographic variables were noted for the selected plots.The stand-structural variables were: canopy cover, shrub cover, basal area, and depth of the litterfall layer.Canopy cover, i.e., the proportion of the forest floor covered by the vertical projection of the tree crowns, is suitable for describing stand-level microclimate and can easily be estimated, thus, most forest inventories record it routinely [32,33].Canopy and shrub covers are expected to be correlated in forest stands [34].The soil and physiographic variables were: texture, content of organic matter, rockiness, slope, altitude, and aspect.The 3SNFI measured cover in percentage, basal area in m 2 ha −1 , litterfall depth in classes of 0.5 cm, and distinguishes among three classes of soil texture (1-sandy, 2-loam and 3-clay), three classes of organic matter content (1-high, 2-moderate and 3-low), and five classes of rockiness (from 1 to 5: 0%, 1%-25%, 26%-50%, 51%-75%, and >75%, respectively).The type of silvicultural treatment observed within the 25 m radius and the presence of damage caused by wild or domestic ungulates were also noted.
Mean values of temperature, precipitation, and potential evapotranspiration [35] of June, July, and August were used to characterize the climatic conditions at each plot during the summer, i.e., the most critical season for P. nigra seedlings survival [9].Additionally, we calculated drought length, as the number of months during which the temperature curve (Ti) is above the precipitation curve (Pi) in the climograph (2Ti > Pi), and an Aridity Index as the quotient between precipitation and potential evapotranspiration during the mentioned months [36].The climatic data were obtained for each inventory plot by an interpolation of the information recorded by the Spanish network of meteorological stations from 1971 to 2000, with 1 km spatial resolution [37].We also obtained radiation data from [38], by uniformly distributing ten sampling points within each region.

Statistical Analysis
Both seedling and sapling abundance at each plot were used to address question (1), i.e., what is the scale of P. nigra regeneration failure, as seedlings and saplings will eventually contribute to maintain stocking in managed stands.However, we only considered data of P. nigra seedling abundance in order to address questions ( 2) and ( 3) as, working with the smallest and presumably youngest individuals, it could be assumed that the stand-structural variables recorded by the 3SNFI surveyors were equal or rather similar to those experienced by seedlings at the time of establishment (P.nigra surpasses 30 cm height in about six years [39]).Therefore, the stand-structural characteristics observed in plots where seedlings are present should be suitable descriptors of the P. nigra regeneration niche.
Data of seedling abundance from the 3226 selected plots were used to address questions ( 2) and (3) as follows.In a first step, we analyzed how seedling abundance varies across levels of canopy cover, using histograms and descriptive statistics.Results suggested that seedling establishment tends to fail in open areas, whilst it is more likely to occur in sites with moderate cover, and decrease in places of high cover (Figure 2).Consequently, we built 2 × 3 contingency tables to test whether regeneration success depends on the level of cover across the study regions: CAT, MAE, IBE, and BET.Plots were cross-classified according to the variables seedlings present / total absence of regeneration, and three levels of canopy cover: low (canopy cover < 55%), moderate (canopy cover: 55%-80%), and high (canopy cover > 80%).These levels were quantitatively established after examining the lower and upper quartiles from the distribution of canopy cover in plots where seedlings were present (Figure 2).In a second step, we generated random numbers to select 170 plots for every combination of region and regeneration level; 170 being the minimum equal number of plots in any of the eight region-regeneration combinations (four regions: BET, IBE, MAE and CAT, and two regeneration levels: seedlings present and total absence of regeneration).This made a randomly selected sample of 1360 plots that was analyzed by Spearman correlation coefficients, Generalized Linear Models (GLM) and multivariate analyses.This selection process was done, as the design needed to be balanced in order to carry out a multivariate analysis of variance (see later).It could be argued that this way of imposing balance dropped part of the available data, but the method is statistically valid [40], and, still, we could analyze data from 136,000 ha of forest.
The relationship between seedling abundance and environmental variables was analyzed by a GLM, with a negative binomial distribution for the error term and log as the link function.A negative binomial distribution was preferred over a Poisson distribution as the response variable showed over-dispersion [41].The predictor variables included in the full model were the stand-structural, soil, and physiographic variables listed earlier, plus the Aridity Index which integrated the observed variability in meteorological data (see later).Canopy cover and shrub cover were included as a linear and a second-order polynomial term to select the best transformation of these two explicative variables to account for non-linearity.The response variable was seedling abundance transformed into a continuous scale by adopting the lowest value of each class range as a conservative approach, i.e., 0, 1, 5, and 16 seedlings per plot for the 3SNFI classes absent, low, medium, and high, respectively [42].We did not detect problems of multicollinearity in the fitted model (Variance Inflation Factor ≤ 1.93 for all the predictor variables).A stepwise procedure was run to select the most significant variables using the Akaike Information Criteria (stepAIC library MASS in R [43]).
Additionally, the effect of woody canopies on the presence or absence of P. nigra seedlings across the four study regions was analyzed by a multivariate approach.This option was chosen, because stand-structural variables can act concomitantly.For instance, canopy cover may facilitate seedling establishment, but accumulates litterfall on the ground, which, in turn, might lessen early seedling survival [9,34].We used a permutational multivariate analysis of variance (PERMANOVA, hereafter).PERMANOVA gives a partitioning of multivariate variation, defined by a distance measure, according to individual factors in any fully balanced multi-way analysis of variance design, with tests done by permutation.No explicit assumptions regarding the distribution of the original variables need to be met [44,45].In addition to PERMANOVA, a canonical discriminant analysis of principal coordinates (CAP) was conducted to test the hypothesis of no differences among the eight combinations of region and regeneration level, and represent them in a multidimensional space graphic.Specifically, the first two axes of the CAP were used to plot an unconstrained MDS (metric multidimensional scaling) of the data.The axis scores were correlated with the original variables using the Spearman correlation coefficient.
Finally, among-region differences for each study variable in plots with seedlings were tested by one-way analyses of variance (ANOVA), or the Kruskal-Wallis test when variables were in a nominal scale.The assumption of homocedasticity in ANOVA was not met for all the variables, even after data transformation.Nevertheless, ANOVAs were still carried out on raw data because they are robust to departures from the assumptions, if data are balanced (i.e., sizes of samples are all made the same) and samples are large [46].The analysis were performed with the computer programs PERMANOVA and CAP [47,48], using 9999 permutations and the Gower dissimilarity as distance measure, because the dataset contained a mixture of continuous and categorical variables [49].GLM and additional analysis were conducted with R 2.11.1 [43].Throughout the paper, values are means ± standard errors.

Climate and Soil Gradients
Results confirmed that the selected four study regions lay along a latitudinal gradient of summer aridity.Radiation did not differ among regions (F 3,36 = 1.25; p = 0.306), but precipitation, mean temperature, potential evapotranspiration and the Aridity Index differed amongst regions for the months of June, July, and August (F 3,1356 ≥ 110; p < 0.0001 for the four variables).Mean temperature and potential evapotranspiration increased, north to south, although the proportion of variability in these variables accounted for by the latitudinal gradient was low (R 2 ≤ 0.24 for both variables), probably, due to an effect of altitude, which increases towards the south (Table 1).However, the latitudinal gradients for the Aridity Index and for the summer precipitation were clearer (R 2 = 0.71 and 0.81, respectively).This north to south gradient of increasing aridity was amplified by the duration of summer drought: 1.3, 2.5, 2.9, and 3.1 months from the northernmost study region to the southernmost one (Figure 1).As a consequence, it is possible to state that P. nigra seedlings are more likely to experience summer drought in localities from the BET region, and less likely in localities from the CAT region.Soil texture was similar in all the study regions, but BET region exhibited the lowest organic matter content (F 1,1358 = 16.59;p < 0.0001), and the highest rockiness (F 1,1358 = 124.04;p < 0.0001).Indeed, rockiness increased gradually towards the south along the latitudinal gradient (β = −0.025;F 1,1358 = 271.6;p < 0.001).

Scale of Regeneration Failure
P. nigra regeneration, either in the form of seedlings or saplings, was present at 52.9% of all study plots, with saplings being present in a higher proportion of plots in the four study regions (Figure 3).The percentage of plots lacking regeneration increased along the rockiness and aridity gradients (30%, 51%, 46% and 58% for the regions CAT, MAE, IBE, and BET, respectively).Multiple logistic regression showed that summer drought, expressed as the Aridity Index, and rockiness were negatively related to the presence of seedlings (df = 1; Wald ≥ 7.01; p ≤ 0.008).Figure 3 shows both a continuous and gradual decline of frequency from plots lacking regeneration towards those with a high density of seedlings or saplings, and a clear preponderance of plots lacking regeneration or belonging to the low class of regeneration abundance.A weak positive correlation was found between seedling and sapling abundances in the CAT, IBE, and BET regions (Spearman's rho = 0.21-0.30;p < 0.003), and moderate in the MAE region (Spearman's rho = 0.47; p < 0.0001), suggesting the existence of a great spatial-variability in the abundance-life stage combinations of regeneration, and that, generally, high densities of seedlings and saplings did not coincide in the same plots.As result, it could be reckoned that regeneration density is at best 509 trees per hectare (four trees/plot) in most forest-stands of the study area.
Table 1.Mean values and standard errors of the abiotic and stand-structural variables measured in 170 plots with presence of seedlings from each study region.Aridity Index = summer precipitation/potential evapotranspiration during summer.Lower values of the index represent drier conditions.The period June-July-August was considered as summer.Different letters between brackets indicate significant differences in a Tukey HDS test (differences were considered significant at p < 0.005 after Bonferroni corrections for 10 variables).*** Significant differences (p < 0.001) in a Kruskal-Wallis test.

Biotic and Abiotic Factors Determining Regeneration
Spearman rank correlation analysis revealed significant univariate relationships between regeneration abundance and several stand-structural, physiographic and climatic variables (Table 2).Table 3 summarizes the frequencies by which seedlings are present or absent across three levels of canopy cover.The observed differences were statistically significant (χ 2 > 8.77; p < 0.004 for all the study regions).Thus, seedlings of P. nigra were more frequent than expected by random in the plots of moderate and high canopy cover and less frequent than expected in the plots of low canopy cover, with the exception of the CAT region, where seedlings were more frequent than expected by random at plots of moderate canopy cover only.Canopy cover explained most of the deviance in the GLM (Table 4).Basal area and shrub cover were also significant explanatory factors of P. nigra regeneration, although the deviance explained by canopy cover was four times larger than the deviance explained by shrub cover (canopy cover and shrub cover were both structural variables that might facilitate seedling establishment).Rockiness and Aridity Index were the only two abiotic factors that influenced P. nigra regeneration.In terms of explained deviance, the contribution of rockiness was the second most important.Coefficient estimates indicated a positive effect of canopy cover, shrub cover and basal area on regeneration, while the effects of rockiness and Aridity Index were negative.For canopy cover and shrub cover a model including the second-order polynomial transformations was selected, indicating a non-linear response of P. nigra recruitment to these two structural variables.The parameter of the second-order term was negative for both canopy cover and shrub cover, suggesting that high values of cover exert a negative effect on seedling abundance, and that seedling abundance might be enhanced at intermediate values of cover (Figure 4).This interpretation was further supported by negatively skewed frequency distributions of canopy cover in plots with presence of seedlings (Figure 5).

Table 4.
Best model (Generalized Linear Model with a log link and a negative-binomial distribution) of Pinus nigra subsp.salzmannii seedling abundance estimated by a stepwise procedure.The selected transformations poly (2) are polynomial second-order transformations, degree of freedom (df), residual degree of freedom (Resid.df), residual deviance (Resid.Dev.) and probabilities χ 2 tests of the effect of the variable are given.Model deviance = 0.064.The sampled plots included a list of 85 shrub species.Among them, the highest percentage of shrub cover corresponded to Mediterranean pioneer plants: genus Thymus, Rosmarinus, and Lavandula, followed by legumes: genus Genista, Erinacea, and Dorycnium, and spiny shrubs: genus Rubus, and Rosa.Respectively, these three groups represented 47%, 23%, and 9% of the overall percentage of shrub cover.An additional GLM, carried out to study the relationship between seedling abundance and cover of the mentioned shrub genuses, showed that facilitation effects are species-specific, being negative for Thymus (df = 1; deviance = 8.99; P(χ 2 ) = 0.002), positive for Lavandula (df = 1; deviance = 5.36; P(χ 2 ) = 0.021) and Rosa (df = 1; deviance = 5.82; P(χ 2 ) = 0.015), and not significant for the other genuses (df = 1; deviances ≤ 2.45; P(χ 2 ) ≥ 0.118).Significant interactions were detected between the region of origin and canopy cover (df = 4; deviance = 38.26;P(χ 2 ) < 0.0001) and shrub cover (df = 4; deviance = 17.89;P(χ 2 ) < 0.01).Seedling abundance was positively affected by canopy cover in the four study regions (Z values ≥ 4.95; P(χ 2 ) < 0.0001).However, the positive effect of shrub cover was limited to the BET region (Z value = 4.06; <0.0001), being not significant for the other three regions (Z values ≤ 1.85; P(χ 2 ) ≥ 0.065).A closer examination of data revealed that legumes were as frequent as pioneer plants in the BET area, and that cover of Erinacea was positively correlated to seedling abundance in this region (Spearman's rho = 0.212; p < 0.001).
The silvicultural treatment applied in 26% of the study plots was selective felling, whilst non-silvicultural treatment was observed in 72% of the plots.Nevertheless, selective felling must be widely implemented in P. nigra forests, as most plots were classified as belonging to multicohort stands by SNFI surveyors (72%, 61%, 79%, and 66% for the regions CAT, MAE, IBE, and BET, respectively).No damage caused by wild or domestic ungulates was observed in the selected plots.

Regional Patterns of Stand-Structural Effects on Regeneration
PERMANOVA detected regional differences in the structural variables considered simultaneously (F = 44.29,p < 0.0001).The structural characteristics of the plots occupied by seedlings differed from those with no regeneration (F = 12.08, p < 0.001), and the interaction region × regeneration density was significant as well (F = 9.26, p < 0.0001).CAP analysis run to test the null hypothesis of no differences among the groups resulting from the region x regeneration interaction further confirmed this result (trace statistic = 0.317, first squared canonical correlation = 0.154; p < 0.0001 for both tests).The first two axes from CAP explained over 87% of the variation observed in the data.The first axis of this ordination showed a strong negative correlation with canopy cover, basal area, and litterfall depth, while the second axis showed a strong negative correlation with shrub cover.These results indicated the existence of gradients along axes 1 and 2, where sites with higher values of canopy and shrub cover were placed at the bottom left side of the graph.In this respect, groups with presence of seedlings tended to dispose themselves around intermediate values of cover and related structural variables, whilst groups representing sites with absence of regeneration exhibited a more disperse pattern (Figure 4).Within regions, seedlings were found to establish in sites with higher values of canopy cover (first CAP axis) in the regions located at mid position along the aridity gradient (MAE and IBE).On the contrary, seedling presence was higher at lower values of canopy cover in the less dry region (CAT), and no differences with regards the effect of canopy cover on the presence or absence of seedlings was found in the driest region (BET).Mean values of canopy cover, shrub cover, basal area, and litterfall depth differed among regions (F 3,1356 ≥ 18.58, p ≤ 0.0001 for the three variables).Thus, CAT region exhibited the highest mean values of canopy cover and litterfall depth, IBE region the highest value of shrub cover, and CAT and BET regions the highest values of basal area (Table 1).

Scale of Regeneration Failure
Our study revealed that regeneration failure of P. nigra might be a common phenomenon across the species distribution area.Nearly half of the plots lacked seedlings or saplings, although the presence of juvenile trees would generally be expected within the observed multi-cohort structure of P. nigra stands.Obviously, the reported patterns of recruitment were influenced by the sample effort applied, just 78 m 2 within a larger plot of 1963 m 2 .Increasing the surveyed area of regeneration in each plot may have provided different recruitment patterns.However, both the area analyzed and the sample size were large enough to accomplish the aims of the study.Thus, results showed a clear decline of regeneration success along the rockiness and summer drought gradients.Seedlings and/or saplings were present in 70% of the plots from the northern region (CAT), but the percentage decreased steadily towards the driest southern edge.Only 42% of the plots from the BET region attained new trees.[42] found significant effects of rockiness on P. nigra recruitment, but, to our knowledge, no other work has previously considered rockiness as an important variable affecting recruitment success in Iberian pine forests.Here, we have found that regeneration is less frequent in plots from the drier and rockier southern edge of P. nigra natural forests, although data did not let us clearly disentangle the contribution of each variable.Other studies have mentioned precipitation as a factor positively correlated with significant increases in recruitment success [42,50].The regeneration pattern observed at a regional scale in this study is then consistent with results from small-scale studies, which identified drought as the main cause of seedling mortality in rock-free P. nigra stands [6,9,27].Consequently, summer drought might be the major limiting factor of recruitment throughout the species distribution area.This impression was further sustained by the facts that, first, overgrazing does not seem to be responsible for the lack of P. nigra regeneration, although herbivore pressure is characteristically high in Mediterranean forests [8,50].Second, although data did not show significant evidences of recent silvicultural treatments (a disturbance that usually triggers regeneration in managed stands), values of stand basal area were within those more suitable for the establishment of new trees [51,52].
Finally, it must be noticed that, where seedlings and saplings are present, current density may be insufficient to restock uneven-aged P. nigra stands.Specifically, the reverse-J-shaped diameter frequency distributions proposed for the species management would need as many as 1600 seedlings per hectare to sustain harvest operations [53].Likewise, at least 2000 seedlings per hectare are needed to restock even-aged stands [51], but we reckoned that current regeneration density is at best 509 trees per hectare.Concerns about the long-term persistence of P. nigra forests under current climate change seem to be justified [26].

Biotic and Abiotic Factors Determining Regeneration
While results indicated that summer drought affects the presence or absence of P. nigra juveniles across the species distribution area, summer drought did not seem to produce any effect upon the density of regeneration where it occurs, as indicated by a non-significant Spearman's correlation coefficient between seedling density and the Aridity Index (Table 2).GLM further confirmed that the density of regeneration is, if at all, only slightly dependent on summer aridity (Table 4).Accordingly, plots with a high density of seedlings and/or saplings were not more frequent at the northern sites than in the southern ones (Figure 3).
Low or non-significant Spearman's r-values, and low explained variances showed that occurrence of regeneration scarcely varies over topographic, and edaphic site factors in the P. nigra distribution area, with the exception of rockiness (Tables 2 and 4).The presence of rocks was associated to poor-sandy soils, steep slopes, and open canopies (data not shown).Therefore, and beyond the fact that rocks diminish the surface available for seed germination, the negative correlation found between seedling density and rockiness could be a consequence of water stress due to a combination of both inferior edaphic conditions, and high evaporation under low canopy cover [8].Canopy cover consequently appeared as the main ecological factor determining seedling establishment in P. nigra stands.

Patterns of Facilitation across P. nigra Distribution Area
Our results indicated that P. nigra juvenile trees can adapt to shady conditions, most likely in order to improve their water balance.Seedlings of this species exhibit no photo-inhibition and a high photosynthetic capacity in full sunlight [54], but recruitment proved to be positively associated with tree cover under natural conditions and across the species distribution area, with the highest frequency of seedling establishment being met at mid canopy cover (between 55% and 80%).Within each study region, the frequency distributions of canopy cover in plots with regeneration exhibited lower coefficients of variation and more negative skew than the plots with no regeneration (Figure 5).Moreover, the environment described by a combination of four stand-structural variables involved in P. nigra recruitment was quite similar for plots with seedling presence across the species distribution area (Figure 4).These findings suggest that cover from the tree canopy is a large-scale habitat condition of the P. nigra regeneration niche, similar to climate or soil conditions.The observation that young P. nigra trees benefit from the special abiotic conditions beneath tree canopies suggests the existence of facilitative interactions between seedlings and canopy trees [6].Indeed, Spanish foresters have been largely aware of the shade tolerance of P. nigra seedlings, and of the importance of tree cover for P. nigra regeneration [51,53], but this study offers insights into the nature of such facilitative interaction.Percentage cover in plots with regeneration did not increase towards the South as we had previously hypothesized based on predictions from the stress-gradient hypothesis.On the contrary, evidence from this study suggests that the interacting effect of light availability and summer drought on P. nigra survival at the seedling stage is nonlinear (humped-back shaped [24]) at both regional and local scales.Overall a humped-back shaped pattern of facilitation was indicated by the inclusion of canopy cover as polynomial of order 2 in the best GLM (Table 4), and by the negative skewness (−0.72 ± 0.07) calculated for the frequency distribution of canopy cover in plots with seedlings present (Figure 2).Within regions, results similarly showed more recruitment at intermediate values of tree cover than at both ends of the canopy gradient (Table 3 and Figure 5).
Most studies about the recruitment dynamics of Mediterranean tree species have identified drought as the major cause of seedling mortality [1][2][3][4], including studies on P. nigra regeneration [6,9].High irradiance increases evaporative demands and further reduces soil-water availability, therefore, there are obvious reasons to explain the increased seedling mortality observed in open areas [5][6][7][8].In poor light environments, the eventual death of P. nigra seedlings because of desiccation can be justified by at least two reasons.First, soil water content may be lower than expected under tree canopies because of rainfall interception and, second, a shift from facilitation to competition for water might occur above a given threshold of canopy cover (~80% in this study) [55].In both cases, seedlings would die in dense stands, because water absorption and light capture cannot be maximized simultaneously, as predicted by the trade-off hypothesis [21] and demonstrated in greenhouse experiments by [54,56] Both reasons could explain why plots with seedling presence from the wetter northern CAT region attained the highest percentages of canopy cover.
However, the trade-off hypothesis [21] and the humped-back shape model [24] proposed a switch from facilitation to competition at both ends of an aridity gradient.This switch was observed in this study in the wettest region but not in the driest one.Specifically, results suggest the existence of competition between seedlings and grown trees in the CAT region: seedling establishment was not more frequent than expected under the more dense tree canopies (Table 3), and regeneration was absent in plots that attained the highest values of canopy cover (Figure 4).On the contrary, sites with and without P. nigra regeneration attained similar values of canopy cover in the BET region, which may suggest a collapse of facilitation in this area with harsher environmental conditions for seedling establishment [25,57].Indeed, a clear facilitative effect of tree canopy on regeneration was detected in the regions located at mid position on the aridity gradient (MAE and IBE), as shown by the distances between groups in the unconstrained MDS of these two regions, representing sites where seedlings are present or absent (Figure 4).
The facilitative effect of shrubs on P. nigra regeneration was not generally extended or strong at the regional scale.Yet, some studies on pines have reported a nurse effect of shrubs.For instance, [6] observed facilitative effects of Juniperus communis on P. nigra regeneration and two sowing experiments carried out in locations within the P. nigra distribution area revealed facilitative effects from shrubs on the seedling survival of Pinus sylvestris, a pine species that also inhabits Mediterranean mountains [4,10].Nevertheless, shrub-tree seedling facilitative effects tend to be species-specific [58], and our results are congruent with the poor nurse quality quoted for the most abundant shrub species in the whole study area [59].Thus, BET was the only region where facilitative effects by shrubs were observed because legume-shrub species were abundant, while the other regions did not benefit from the poor nurse quality of the more frequent Mediterranean pioneer species.

Conclusions
This study identifies that P. nigra may have recruited insufficiently during the last decades in order to restock managed stands, this circumstance being more pronounced in the drier southern edge of the species distribution area.The canopy of adult trees facilitates seedling establishment, preferably within an interval of 55%-80% canopy cover.Results suggests that this facilitative interaction is non-linear (humped-back shaped) at the local (within region) and regional (across regions) scales, although this pattern interacts along the aridity gradient showing competition in the wettest site, high facilitation in the mid-dry sites, and low facilitation at the driest site.This pattern resembles predictions from the stress-gradient hypothesis, but with a collapse of facilitation in the most water-stressed environment.We detected these patterns by analysing data along a continuous aridity gradient (GLM), and among regions (contingency tables and multivariate analysis).Evidence suggests that some shrub species may facilitate recruitment in the harsher areas.Collectively, our results reduce the possibilities of adapting forests management to drought by the application of alternative silvicultural prescriptions involving canopy cover.Foresters working with Mediterranean trees at the southern limit of their distribution area will need to know as much about shrub ecology as they already know about the ecology of commercial tree species.

Figure 1 .
Figure 1.Geographical distribution of Pinus nigra ssp.salzmannii in Spain with indication of the four study regions: Prepirineo Catalá n (CAT) in yellow, Alto Maestrazago-Montes Universales (MAE) in red, Sistema Ibé rico Meridional (IBE) in blue, and Cordilleras Bé ticas (BET) in green.The inset indicates the distribution area of Pinus nigra sensu lato.

Figure 2 .
Figure 2. Frequency distribution of canopy cover in plots with seedlings present (n = 913).Arrows indicate the position of different percentiles.

Figure 4 .
Figure 4. Canonical discriminant analysis of principal coordinates of the stand-structural variables characterizing the sites where P. nigra subsp.salzmannii seedlings are present (pre) or absent (abs) in the four study regions (CAT, MAE, IBE, and BET).Data are means ± standard errors (n = 1360).

Figure 5 .
Figure 5. Frequency distribution and descriptive statistics of canopy cover in sites where seedlings were present (A) or absent (B) within each study region.SD is the standard deviation; CV is the coefficient of variation.

Table 2 .
Spearman rank correlations between biotic and abiotic variables, and seedling density across the four study regions; significant results are noted in bold.