Ecological Diversity within Rear-Edge: A Case Study from Mediterranean Quercus pyrenaica Willd.

: Understanding the ecology of populations located in the rear edge of their distribution is key to assessing the response of the species to changing environmental conditions. Here, we focus on rear-edge populations of Quercus pyrenaica in Sierra Nevada (southern Iberian Peninsula) to analyze their ecological and ﬂoristic diversity. We perform multivariate analyses using high-resolution environmental information and forest inventories to determine how environmental variables differ among oak populations, and to identify population groups based on environmental and ﬂoristic composition. We ﬁnd that water availability is a key variable in explaining the distribution of Q. pyrenaica and the ﬂoristic diversity of their accompanying communities within its rear edge. Three cluster of oak populations were identiﬁed based on environmental variables. We found differences among these clusters regarding plant diversity, but not for forest attributes. A remarkable match between the populations clustering derived from analysis of environmental variables and the ordination of the populations according to species composition was found. The diversity of ecological behaviors for Q. pyrenaica populations in this rear edge are consistent with the high genetic diversity shown by populations of this oak in the Sierra Nevada. The identiﬁcation of differences between oak populations within the rear-edge with respect to environmental variables can aid with planning the forest management and restoration actions, particularly considering the importance of some environmental factors in key ecological aspects.


Introduction
The study of ecological dynamics within the rear-edge populations is considered essential to establish proper management guidelines under current climate uncertainties [1]. Rear-edge populations are often adapted to local environmental conditions at the limit of the species' ecological amplitude, and often show a long-term persistence [2]. Local responses to environmental changes may differ from the species mean response [3][4][5][6], and such differences may either promote or hamper the survival of edge populations under global change [7]. Furthermore, the heterogeneity in the response to climate change observed across ecological and geographical gradients [8][9][10][11], justifies the need to incorporate fine-scale variation of environment variables throughout species ranges to better understand species responses to global change [12,13]. This is particularly important for mountain landscapes, where the topographic complexity may cause a decoupling between the climate and the geographic spaces [14,15].  Table 1.

Quercus pyrenaica Forests
Quercus pyrenaica is a deciduous species extending through southwestern France, the Iberian Peninsula, and northern Morocco [28]. The forests of this species reach their southernmost European limit in Andalusian mountains such as Sierra Nevada, where eight populations have been identified (Figure 1a; Table 1) on the basis of their isolated geographic locations in deep valleys separated by distances considerably longer than the average dispersal distances of the seeds by birds such as the Eurasian jay (Garrulus glandarius) [41,42]. They are distributed on siliceous soils both in the northwestern and southern slopes of the mountain range and are often associated with major river valleys. These oak woodlands represent a rear edge of their distribution [2], containing high levels of intraspecific genetic diversity [43]. Their conservation status for southern Spain is "Vulnerable", and it is expected to suffer from climate change, potentially reducing its suitable habitats in the near future [8,44]. The distribution of Q. pyrenaica forests in Sierra Nevada was delimited using the updated version of the forest map of Sierra Nevada on a 1:10,000 scale [40,45]. Black and white ortophotographies from 2001 (0.5-m of spatial resolution) and false color aerial photographies (Color Infrared) from 2005 (1-m resolution) were used to correct errors by detailed photographic interpretation, resulting in a detailed map of oak forests (Figure 1b). Forest patches with at least 50% tree cover of which 75% cover being Q. pyrenaica were considered oak patches.

Environmental Data
For each oak population, we obtained the values of 30 environmental variables selected to represent different direct and indirect gradients important for plant distribution [46,47]: temperature, water availability, topography, solar radiation, and land-use (Table 2). Observed climate data (1960-2010) from 43 meteorological stations 50 km around Sierra Nevada, compiled by Sierra Nevada Global Change Observatory [21], were used as input to compute high resolution (100 × 100 m pixel-size) climate maps [48] based on the mapping method proposed by Ninyerola et al. [49]. Seasonal and annual maps with the averages of direct solar radiation and insolation time were computing using the GIS GRASS module r.sun [50,51]. From a high-resolution digital elevation model (10-m; Department of the Environment, Regional Government of Andalusia), several topographic variables were derived: elevation, slope, aspect, E-W, and N-S gradients, topographic position (difference in elevation between a cell and surrounding cells within a 1000 m radius) [52]. In addition, topographic wetness index and flow accumulation were computed using the r.terraflow module of GRASS GIS. As a surrogate of anthropogenic influence, we computed the frequency of human infrastructures in a 2000 m radius buffer. Finally, for each environmental variable, we extracted the values for all the 100 m size pixels contained within each oak population ( Figure 2). Using an environmental data matrix, the main environmental gradients that characterize the oak forests at Sierra Nevada were identified using a Principal Component Analysis (PCA). Linear Discriminant Analysis (LDA) was also applied to identify different groups of oak populations. With a matrix of floristic composition, a Nonmetric Multidimensional Scaling (NMDS) ordination were applied to visualize patterns of species composition, interpret them according to the environmental factors, and identify groups of oak populations based on similarities between floristic composition. See Materials and Methods for more details.

Forest Attributes
To characterize oak patches, we selected several stand attributes relating to forest structure, function, and composition from Sierra Nevada Forest Inventory [53] (Table 2). By using this approach, we characterized the plant community both in terms of their species composition, and also regarding their ecological functioning [54,55]. SINFONEVADA forest inventory was carried out during 2004-2005, and it includes an extensive network of plots distributed within the main forest units of the Sierra Nevada mountain range. We selected 32 plots belonging to the deciduous broadleaf forests category. All of them are located within the eight Pyrenean oak populations identified in Sierra Nevada. For each plot (20 × 20 m), all trees with diameter at breast height (dbh) >7.5 cm were tallied by species and dbh. Regeneration, species composition, and abundance were also recorded in two additional subplots (see [53] for a detailed description): a 5-m radius subplot where the seedling abundance of Q. pyrenaica was recorded; and a 10-m radius subplot where the species composition and abundance estimated by the Braun-Blanquet cover-abundance scale were measured [56] (see Table S1).
Forest composition (richness) and plant diversity were used as an indicator for overall forest biodiversity. Plant diversity was measured using the Shannon diversity index [57]. The total regeneration was used as a proxy for forest functioning. Finally, as forest structure indicators, we selected the following attributes: the total-and strata-(i.e., tree, shrub, and herbaceous) canopy cover; canopy cover diversity; tree height, tree density, basal area, and volume of adult tree. Canopy covers were computed as the proportion of plot area covered by the whole forest (total) and the different strata considered (tree, shrub, and herbaceous, respectively). Canopy cover diversity was quantified through the Shannon index for the proportion of plot area covered by different vegetation strata (tree, shrub and herbaceous) according to the following equation: CCd = ∑ n i=1 g i · ln g i , where g i is the proportion of strata i of the total plot area and n is the number of strata [58]. Basal area was calculated as the sum of the basal areas of the adult trees assuming a circular cross-section of the trunk. Volume was calculated as the sum of volume (V = 0.55 × height × diameter 2 ) of all Q. pyrenaica adult trees. Additionally, we also extracted the values of the environmental variables for the centroids of the plots, and we added a species-composition matrix for each of the 32 selected plots.

Statistical Analysis
To identify the main environmental gradients that characterize the oak forests at Sierra Nevada, we performed a principal component analysis (PCA) on the standardized variables ( Figure 2). Over 75% of the correlations (Sperman's r) among variables were significant (p < 0.01). We checked the adequacy of the environmental matrix by applying the Kaiser-Meyer-Olkin test, a measure of sampling adequacy (KMO = 0.7138, value greater than 0.5 is considered adequate [59]. The Kaiser-Guttman rule [60], i.e., axes whose eigenvalues are larger than the average of all eigenvalues; and the criterion that any principal component (PC) accounts for at least 10% of the total variance were used to determine the meaningful PCs to be retained for interpretation [61]. The PCA variables with a correlation to the principal components that was higher than 0.7 were selected to describe the environmental gradients indicated by the principal factors. We applied Linear Discriminant analysis (LDA) to determine the environmental variables that best discriminated among Pyrenean oak patches and to identify different groups of populations [61,62].
Then, environmental variables and forest attributes were tested for differences among population groups previously identified. Normality and homoscedasticity were checked using the Shapiro-Wilk test and Levene's test, respectively. If normality and homoscedasticity assumptions were satisfied, we performed ANOVA analysis followed by the Tukey LSD for testing statistical significance. Otherwise, Kruskal-Wallis ANOVA for nonparametric data were conducted followed by manual pairwise comparison using the Mann-Whitney U-test.
Finally, we used a Non-metric Multidimensional Scaling (NMDS) ordination analysis based on Bray-Curtis dissimilarity distance [63] to: (i) visualize patterns of species compositions, (ii) interpret them with respect to the environmental factors (i.e., relate the variability in species composition to environmental variables), and (iii) identify groups of Pyrenean oak populations based on similarities between floristic composition. NMDS involves the reduction of multidimensional similarity data to a low-dimensional ordination in which relative distance indicates relative similarity (i.e., plots with very similar species composition are close and vice versa) [64]. We compared two and three-dimensional solutions based on Kruskal's stress (as a measure of goodness of fit). We also studied the floristic-environment relationships by fitting linear trends on the ordination yielded by the NMDS. For these linear fittings, squared correlation coefficients and empirical p-values were calculated using random permutations (n = 1000) of the data [65]. Finally, we fitted non-parametrically smoothed surfaces of continuous environmental variables on the NMDS ordination. The smooth surfaces were fitted using generalized additive models (GAM) with thin plate splines, using the coefficient of determination (R 2 ) as goodness-of-fit statistics e.g., [65,66].

Results
PCA of all measured environmental variables that yielded three significant axes explained 62.11% of the total variance ( Table 3). The first PC axis was strong and negatively correlated with radiation and precipitation related variables, and positively with northness gradient and slope (Table 3). Maximum average temperatures showed the strongest negative correlations with the second PC axis. The third PC axis was negatively correlated with minimum average temperatures. The precipitation variables presented weak positive correlation with the third PC axis.
The discriminant analysis yielded three significant functions explaining 97.9% of variance ( Table 3). The ordination plot ( Figure 3) showed a clear separation of oak populations into three clusters: a single-oak-population (CAM) cluster, namely N in the Figure 3; the second cluster (NW) formed by the GEN, MON, DIL, and DUR oak populations; and the southern cluster (S) composed by the southern oak populations CAN, POQ, and TRE. Southern oak populations were separated out from northern populations along the first LDA axis (Figure 3), which showed slight negatively correlation with autumn rainfall. The second and third LDA axes showed weak correlations with all variables ( Table 3).
The three oak clusters showed significant differences for most of the environmental variables analyzed (Table 4. Only winter minimum temperatures (χ 2 = 5.35; p-value = 0.069) and insolation time during summer (χ 2 = 0.306; p-value = 0.306)) was similar among the three oak clusters (Table 4). Post-hoc analysis showed that, for most of the environmental variables, we found pairwise significant differences between all three of the oak clusters (Table 4).
Forest attributes did not significantly differ among the above described oak clusters except for plant diversity and herbaceous canopy cover ( Table 4). The N cluster showed a higher value of Shannon diversity index (2.27 ± 0.17) than the NW cluster (Mann-Whitney U = 22.0; p-value <0.01). For stand attributes relating to forest structure, only the herbaceous canopy cover showed significantly differences (χ 2 = 11.18; p-value = 0.004; Table 4) between N and NW clusters (Mann-Whitney U = 15.0; p-value < 0.01). For all other forest structure attributes, despite there being no significant differences, the N cluster showed the lowest values (Table 4). No significant differences were recorded for regeneration variables.
A three-dimensional solution of the NMDS was chosen because its correlation with the original data was higher than for a two-dimensional solution (Linear fit R 2 = 0.793 vs. 0.713). Additionally, lower Kruskal's stress value was observed for the three-dimensional solution (Stress = 0.159 vs. 0.226). The NMDS ordination of the forest stands according to their floristic composition was significantly correlated with precipitation variables, elevation, and marginally with winter maximum temperatures ( Figure 4; Table 5). The precipitation variables showed highly and negative correlations with NMDS axis 2 ( Table 5). The NMDS axis 1 were negatively correlated with elevation (R 2 = 0.464) and minimum temperatures, and positively correlated with slope and winter maximum temperatures ( Table 5). The NMDS ordinations with fitted vectors and surfaces for significant variables are shown in Figure 5. All of these variables showed a nonlinear significant relationship with the ordination pattern (R 2 values for surfaces were slightly higher than linear R 2 values; Table 5).   Table 1. Numbers in brackets expressed explained variance (%) for each discriminant axis.   Table 5). Each plot was colored according to the three oak-populations' clusters derived from discriminant analysis. Only two dimensions of the NMDS are illustrated for ease of representation.

Ecological Diversity within the Rear-Edge
The rear-edge populations of Quercus pyrenaica located in mountain areas are not ecologically homogeneous, neither for their environmental conditions nor for their plant species composition. In this study, we find separate groups of Q. pyrenaica populations within Sierra Nevada (rear-edge) driven by radiation and rainfall as main discriminant variables (Figure 3). The differences among populations based on environmental variables obtained in our study are in line with differential ecological dynamics reported for Q. pyrenaica forests in the Sierra Nevada by other studies. For instance, primary productivity of these forest measured using remote sensing showed a heterogeneous spatial behavior, with oak woodlands of the southern slopes displaying a greater annual vegetation greenness than those from the northern slopes [11,77,78]. In addition, differences have been found in both seasonal dynamics of greenness [77], and in temporal trends for primary productivity in the last few years related with differential snow-cover trends in contrasting slopes [78,79].
Interestingly, our results also showed differences in species diversity among population groups derived from clustering based on environmental variables. These results are consistent with those provided by Lorite et al. [80], who pointed out that differences observed for the floristic component in the Q. pyrenaica populations of Sierra Nevada are related to the microclimatic conditions. Thus, Lorite et al. [80] found that the oak woodlands located in the northern part of the Sierra Nevada showed greater floristic similarity with those located at the center of the Q. pyrenaica distribution than those located at southern slopes of Sierra Nevada (geographically closer) [80]. The floristic differences between Sierra Nevada oak populations could also be related to the anthropogenic impact suffered by those populations, since the anthropic disturbances can affect the floristic patterns of the woodlands of this species, as it has been documented for oak woodlands in central Spain [81]. Thus, our results shown that the CAM oak population (N-cluster) showed both the highest plant species diversity and richness (Tables 4 and S1), which may be related to a better conservation status, as this population has been less exposed to intense anthropogenic activity [82]. Conversely, the southern oak populations (CAN, POQ and TRE) showed a poorer floristic composition conditioned by both climate and intense land use [83,84].
We found a remarkable match between the population's clustering derived from the analysis of environmental variables ( Figure 3) and the ordination of the populations according to species composition (Figures 4 and 5). These findings suggest a linkage between the heterogeneity of environmental factors and the variability of species composition for these woodlands. The diversity of ecological conditions for Q. pyrenaica populations in this rear edge are in line with the high levels of genetic diversity shown by populations of this oak in the Sierra Nevada [43,85,86]. The oak woodland of Sierra Nevada has shown higher values of both genetic diversity and allelic richness than those populations located in Central Spain [85,86]. For Sierra Nevada oak populations, a great genetic differentiation among populations has been reported. [86]. Specifically, Valbuena-Carabaña and Gil [86] found high values of population genetic differentiation between oak stand located in the El Camarate site (CAM population), on the northern slope of Sierra Nevada, and other stands located in the Cáñar site (CAN population), located on the southern slope. The climatic and topographical heterogeneity that exists in the Sierra Nevada offers a great diversity of microhabitats, which has allowed this mountain range to act as a refuge for different species [87][88][89], including for deciduous Quercus species during the last glacial period [90][91][92]. In fact, there is fossil and genetic evidence for different Quercus species that strongly suggests that they survived only in southerly refugia during the last glacial maximum [90,[93][94][95]. The persistence in a refugium suggests a combination of a moderately suitable local environment buffering against the regional climate, and a relative tolerance to climate change, by either pronounced phenotypic plasticity, and/or adaptive capacity [96]. This could be very well the case of Q. pyrenaica, a species harboring a high genetic diversity [43], located in a mountain region with a complex topography that could protect local populations against rapid climate shifts and allow species to persist despite regionally unfavorable environments.

The Importance of Summer Rainfall at the Micro-Habitat Level
The distribution of Q. pyrenaica is known to be conditioned by summer drought period with a minimum of 100-150 mm of summer rainfall [97,98]. del Río et al. [99] in a bioclimatic analysis for this species revealed the importance of rainfall and ombrothermic indexes in the separation of temperate and Mediterranean forests [99]. On a more detailed scale, the distribution for this oak is driven by a complex gradient related with temperature, rainfall, and radiation [36,100]. Our study unveils a separation in the environmental space between oak populations at the rear edge related with the spatial pattern of precipitation for this mountain region [101]. Thus, summer and annual rainfall are among the most important factors in explaining the distribution of Q. pyrenaica forests in Sierra Nevada ( Table 3). The northern and northwestern populations of Q. pyrenaica in the Sierra Nevada are located in valley bottoms with northern orientation, where the relative humidity is greater as a result of a lower solar radiation. On the other hand, the populations of the southern slopes of Sierra Nevada get an extra supply of water from moist air from the neighboring Alborán sea [102]. The differences in water availability among oak populations could affect several ecological processes such as tree-growth [5,11], seedling germination and survival [103][104][105], and the regeneration of the species [106], mainly due to the key role of water availability in the microsites facilitating the germination and establishment of seedlings.

Implications for Forecasting and Modeling
The factors controlling species distributions may vary depending on the scale of observation. In large scale areas, the distribution of a species is likely to be controlled by climatic regulators [107], whereas, on local scales, factors related to biological interactions play a relevant role in shaping species distributions [108,109]. At the site level, we found that moisture availability is the environmental factor that better separates the studied oak populations into clearly differentiated clusters. The identification of different population groups based on environmental variables at fine-scale is important when modeling the distribution and forecasting the impact of global change on the species. Our results suggest that incorporating the local adaptations of individual populations into predictive models might help avoid misrepresenting the potential range shift of species under changing climate conditions [7]. This is particularly important for species with rear edges located in mountain ranges, since these areas provide a broad diversity of microhabitats due to climatic and topographical heterogeneity [87]. For instance, some recent works have performed high-resolution models of the distribution of relict trees in Mediterranean southern mountains (e.g., Abies pinsapo, Pinus sylvetris and P. nigra) providing useful information for forest management actions [110].

Conclusions: Biodiversity from the Genetics to the Landscape
We identified several groups of oak populations within the rear-edge of the Q. pyrenaica forest mainly due to microhabitat conditions. The different clusters of oak populations are supported both by discriminant analysis of environmental variables and by ordination analysis based on the floristic composition in the target populations. Our results show that the diversity in the ecological conditions within these populations results from both the environmental heterogeneity created by the slopes and the contrasting exposures of the valleys they inhabit, and also the anthropic use of these ecosystems e.g., [11,111]. The confluence of these factors generates a multitude of environmental conditions on a fine scale, which are reflected in the distribution, composition, and functioning of the Quercus pyrenaica forests. Quercus pyrenaica woodlands are highly diverse at all organization levels, from a genetic perspective, i.e., high levels of genetic differentiation within species [43] and differences between populations [86]-to ecosystem-functioning level, i.e., diversity in terms of primary production and growth [78,79], and diversity of resilience to disturbances e.g., [11]. Such ecological heterogeneity is also made evident by the accompanying plant communities, which are very different depending on the oak population considered, such differences being correlated with the differences in environmental conditions among populations.
Mountains such as Sierra Nevada do not only act an elevation gradients along which plant communities are distributed and replaced, in fact, they constitute an ecological mosaic in which other factors besides elevation, e.g., the exposure and the history of human management, create a broad range of responses from the oak woodlands and its very diverse associated vegetation, from genetics to landscape. Understanding the differences that exist between oak populations within the rear edge with respect to environmental variables help us to plan both the forest management and restoration actions, especially taking into account the importance of some environmental factors in key ecological aspects e.g., regeneration and growth [11,104]. Our results also show the importance of the rear-edge mountain areas as a refuge for within-species diversity, and the role of species' southern ranges as hotspots of within-species diversity [112,113]. All of this knowledge will be important to prioritize the conservation measures, and to design adaptive management actions targeting these populations, in order to maintain their ecological processes and biodiversity.