Niche Characterization of Shrub Functional Groups along an Atlantic-Mediterranean Gradient

: The identiﬁcation of the factors controlling the understory species distribution and abundance is essential to understand the ecology and dynamics of natural forests and their management response. We assess the relationships between environmental gradients and shrub functional groups distribution patterns and niche characteristics in a transitional area between the Eurosiberian and Mediterranean biogeographic regions in Northern Spain. Here, 772 plots from the 3rd Spanish National Forest Inventory were used. Shrub functional groups respond to the same complex environmental gradients as trees, i.e., the north-south climatic gradient and a slope gradient. Unimodal response curves of shrub functional groups and families dominate along both gradients, providing evidence of successful functional turnover. Similar to tree species, the niche location of functionally related shrubs is close. Functional groups occupying environments with sharp contrast or transitional environments have the broadest niches, whereas those specialized functional groups occupying localized habitats showed the narrowest niches. The knowledge of shrub species distributions and niche characteristics along complex environmental gradients will improve our ability to discuss potential conservation management goals or threats due to land-use changes and future climate change.


Introduction
Shrub species living in the forest understory are important elements of ecosystem structure and function, providing habitat and forage for wildlife [1][2][3][4] and contributing significantly to plant diversity [5][6][7]. Being particularly important in areas where shrub understory is significantly more diverse than forest overstory, such as the Palencia province (Northern Spain) [8,9]. It is widely recognized that shrubs facilitate the growth of protected species under their canopy [10], besides the recruitment of tree seedlings mostly from Quercus L. (oaks) [11][12][13][14] and Pinus L. (pines) species [15][16][17] in Mediterranean environments. These facilitative effects are fundamental to regenerate forest ecosystems in Mediterranean regions [18,19]. Additionally, shrubs play a fundamental role as ecosystem engineers to increase heterogeneity in micro-environmental conditions favoring late-successional species establishment [12,[20][21][22]. Conversely, it has been shown that certain shrubs favor exotic species invasion (e.g., in the Northern Californian coastal dunes Ericameria ericoides (Less.) Jeps. enhances the colonization of the invasive Bromus diandrus Roth [19]), and even some shrub species become invasive in the lack of appropriate management practices (e.g., Cytisus scoparius (L.) Link in France and Australia [23,24]). Thus, the description of environmental factors shaping the distribution and abundance of understory shrub species has fundamental implications for forest diversity conservation [25,26] and forest management [27,28], particularly in the current scenario of climate change [29][30][31].
Despite the fact that shrubs are key components of forest ecosystems [32,33], there have been few efforts to model shrubs abundance or distribution [5,34], and very few models of simulation of forest dynamics incorporate the response of understory vegetation layers [26]. Most likely, the understory responses are more complex than overstory responses. For instance, microclimatic extremes are ameliorated by moderate overstory cover, which would otherwise be stressful for ground-layer plants, in contrast, high overstory cover suppresses understory species irrespective of environmental constraints [5]. Thus, it might be assumed that understory vegetation has a unimodal response to forest biotic variables, such as overstory cover [5], along with environmental gradients [35,36]. Nevertheless, as reported for tropical forests [37,38], understory shrub species might not respond to the same complex environmental gradients as tree species. At the same time, shrub community responses are often highly complex [39]. Therefore, to reduce complexity, shrub species are usually grouped according to the traits they shared into different functional groups or plant functional types [40], such as Raunkiaer's life-forms (chamerophyte, phanerophyte), dispersal mode (anemochory, zoochory . . . ) or regeneration method (germinator, resprouter). In any case, a deeper understanding of shrub functional responses to complex environmental gradients is needed.
In this sense, we assessed the relationships between shrub functional group distribution patterns and niche characteristics and two main environmental gradients (coenoclines) previously identified [41] in a transitional area between the Eurosiberian and Mediterranean biogeographic regions in Northern Spain. We characterized the realized niche of shrub functional groups along both coenoclines, considering that the species realized niche is conceived as a multidimensional response surface and is composed by the fundamental requirements of the species and its function, i.e., interactions with other species through competitive exclusion and facilitation [42]. Here, considering that functional groups are important components of ecosystem function [43], we used a shrub-functional-groups approach. Our specific objectives were: (1) to assess whether understory shrub functional groups responded to the same regional complex environmental gradients as tree species (i.e., altitude and temperature, and steepness); (2) to determine whether shrub functional groups exhibited a unimodal response to those environmental gradients; (3) to quantify the optima and niche widths of shrub functional groups across the main gradients. The description of shrub functional groups distribution and niche features along complex environmental gradients will help to improve our knowledge to discuss potential conservation management threats and goals because of land uses and climate changes [15].

Study Area and Data Source
The environmental gradient under study runs from the north to the south of the Palencia province (Northern Spain; 43 • 04 N and 41 • 46 N latitude, and 3 • 53 W and 5 • 02 W longitude; Figure 1A). The gradient is 140 km in length and comprises a high variety of environmental conditions and landscapes [9,44], as a result of the confluence of the Atlantic and Mediterranean biogeographic regions and two geomorphological units (the Cantabrian Range and the Castilian plateau), which harbor great vegetation diversity [8]. The climate changes from Alpine to the Mediterranean from north to south [45], which is related to topography ( Figure 1B), mainly determined by the presence of the Cantabrian Range in the north (altitude up to 2540 m), and the Castilian plateau (mean altitude ca. 800 m) in the Centre and southern parts of the gradient. The temperature increases markedly from north to south followed by a notorious increase in precipitation (higher xericity) that leads to a rise in the continental nature of the climatic conditions [46]. Thus, the climatic and environmental gradient throughout the Palencia province is of particular interest when describing and characterizing the forest vegetation compositional changes from spatial perspectives along broad heterogeneous environmental gradients. 800 m) in the Centre and southern parts of the gradient. The temperature increases markedly from north to south followed by a notorious increase in precipitation (higher xericity) that leads to a rise in the continental nature of the climatic conditions [46]. Thus, the climatic and environmental gradient throughout the Palencia province is of particular interest when describing and characterizing the forest vegetation compositional changes from spatial perspectives along broad heterogeneous environmental gradients. Shrub species data from 772 permanent field plots ( Figure 1A) from the 3rd Spanish National Forest Inventory (3SNFI; 1997-2007) were used. The SNFI measure circular plots of variable radius (5, 10, 15, and 25 m) systematically distributed and located on the intersection nodes of a 1 × 1 km 2 Universal Transverse Mercator grid; only plots located inside forest areas are measured (see [41,46] for further details on plots selection, characteristics, and distribution). Trees with a diameter at breast height of 75, 125, 225, and 425 mm are measured, respectively, in each one of the concentric circumferences. The cover (%) of all Shrub species data from 772 permanent field plots ( Figure 1A) from the 3rd Spanish National Forest Inventory (3SNFI; 1997-2007) were used. The SNFI measure circular plots of variable radius (5, 10, 15, and 25 m) systematically distributed and located on the intersection nodes of a 1 × 1 km 2 Universal Transverse Mercator grid; only plots located inside forest areas are measured (see [41,46] for further details on plots selection, characteristics, and distribution). Trees with a diameter at breast height of 75, 125, 225, and 425 mm are measured, respectively, in each one of the concentric circumferences. The cover (%) of all shrub species present in a fixed plot radius of 10 m [46] was inventoried. A total of 86 shrub species were registered in the 772 selected plots. However, only 47 of them were used in subsequent analyses for having sufficient cover and frequency in the plots, a correct taxonomical identification in the field [41], and available information of their functional traits in databases for Southern European Flora. These 47 shrub species were classified in habitat-related functional groups (see Table A1) such as family, geographical distribu- tion (Atlantic, Mediterranean, Atlantic-Mediterranean, endemic of the Iberian Peninsula), biotype (perennial, deciduous), Raunkiaer s life-form (chamerophytes, phanerophytes), dispersal mode (anemochory, zoochory, authocory, barocory, other), and regeneration method (germinator, resprouter).

Coenoclines Characterization
In previous studies, an indirect ordination technique (detrended correspondence analysis, DCA) was applied on the cover (%) matrix of all woody species (trees and shrubs) present in each of the 772 selected plots [41,46]. DCA showed that changes in the forest vegetation compositional changes along the Palencia province were mainly determined by the north-south climatic differences (DCA1, primary coenocline), although a secondary coenocline also showed a turnover of species related to the slope gradient (DCA2, steepness). In particular, Fagus sylvatica L. and Quercus petraea (Matt.) Liebl. dominated deciduous mountain forests that are replaced by Quercus pyrenaica Willd. forests and those, in turn, by Quercus faginea Lam. and Quercus ilex subsp. ballota (Desf.) Samp. Forests, as aridity increases towards the south. Pinus sylvestris L. and Pinus nigra J.F. Arnold dominated coniferous woodlands in mountain areas and detrital moors (middle part of the province; Figure 1B), whereas the southern limestone moors ( Figure 1B) are dominated by plantations of Pinus halepensis Mill., Pinus pinea L. and Cupresus sempervirens L. Interestingly, along the slope gradient there is also a turnover of tree species, parting natural Juniperus spp., Quercus petraea, Quercus pyrenaica, and Fagus sylvatica forests located on steeply sloping sites, and pine plantations (Pinus nigra, Pinus pinaster Aiton) that dominate in flat areas. Therefore, the search for patterns that relate tree species trends with understory shrub functional group trends along these two main coenoclines will give us valuable information for management and conservation.

Data Analyses
First, the main two coenoclines were previously identified (DCA1, DCA2) [41,46]. Therefore, here we modeled the abundance (cover percentage) of shrub functional groups along these coenoclines using HOF models (HOF = Huisman-Olff-Fresco [47]) through the 'eHOF' package (version 1.7 [48]) in the R-language environment (version 4.0.3; R Development Core Team, Vienna, Austria, http://www.r-project.org, accessed on 1 February 2020). HOF models allow describing the species responses, which may result from both environmental conditions and intra-and interspecific interactions [49]. HOFs are a hierarchical set of response models, ranked by their increasing complexity (Model I, no species trend; Model II, increasing or decreasing trend; Model III, increasing or decreasing trend below maximum attainable response; Model IV, symmetrical response curve; Model V, skewed response curve). The selection of the most appropriate model for each of the functional groups analyzed was done using the Akaike Information Criterium (AIC) [50]; smaller AIC values indicate more parsimonious models. Finally, for those functional groups with unimodal responses, the location of niche optima (µ) and niche widths (2t) were derived from the HOF models [49].

Distribution Patterns of Main Shrub Families along Coenoclines
Fabaceae and Caprifoliaceae were the unique families with indeterminate response curves (HOF model I) since both had low and constant cover (<1%) along increasing aridity DCA1 coenocline. Thymelaeaceae and Oleaceae (HOF model II with a decreasing trend) had the greatest cover in the northern mountains with a subsequent reduction towards the south (DCA1 left-end; Figure 2A). Only Asteraceae (HOF model II with an increasing trend) showed an increasing trend as aridity increases towards the south (DCA1 right-end). The remainder taxonomic groups showed unimodal response curves (HOF models IV or V; Figure 2A and Table 1) with optima at different points along the gradient, for example, Aquifoliaceae and Betulaceae with optima in the northern mountains (DCA1 < 2) had the narrowest niche width (2t = 0.5), Ericaceae had intermediate niche width (2t = 2.3) with its optimum in the transitional area between the Atlantic and Mediterranean conditions (DCA1 ≈ 2), Rosaceae and Cistaceae had the broadest niche widths (2t = 6.1 and 4.1, respectively) with optima in the gradient middle (DCA1 ≈ 4), and Labiatae had intermediate niche width (2t = 2.6) with its optimum at the right end of DCA1, i.e., towards the south.
Rosaceae, and Ericaceae accounted for 49% of shrub genera and 44% of shrub species.

Distribution Patterns of Main Shrub Families along Coenoclines
Fabaceae and Caprifoliaceae were the unique families with indeterminate response curves (HOF model I) since both had low and constant cover (<1%) along increasing aridity DCA1 coenocline. Thymelaeaceae and Oleaceae (HOF model II with a decreasing trend) had the greatest cover in the northern mountains with a subsequent reduction towards the south (DCA1 left-end; Figure 2A). Only Asteraceae (HOF model II with an increasing trend) showed an increasing trend as aridity increases towards the south (DCA1 rightend). The remainder taxonomic groups showed unimodal response curves (HOF models IV or V; Figure 2A and Table 1) with optima at different points along the gradient, for example, Aquifoliaceae and Betulaceae with optima in the northern mountains (DCA1 < 2) had the narrowest niche width (2t = 0.5), Ericaceae had intermediate niche width (2t = 2.3) with its optimum in the transitional area between the Atlantic and Mediterranean conditions (DCA1 ≈ 2), Rosaceae and Cistaceae had the broadest niche widths (2t = 6.1 and 4.1, respectively) with optima in the gradient middle (DCA1 ≈ 4), and Labiatae had intermediate niche width (2t = 2.6) with its optimum at the right end of DCA1, i.e., towards the south.  Table 1, with the additional: Thymelaeaceae = Thymel; Oleaceae = Olea.   Table 1, with the additional: Thymelaeaceae = Thymel; Oleaceae = Olea. aceae (HOF model II with an increasing trend) increased their cover towards the steepest sites (DCA2 left-end). The other taxonomic groups showed unimodal response curves (HOF models IV or V; Figure 2B and Table 1) with optima at different points along the gradient: Betulaceae, Aquifoliaceae, and Rosaceae towards the steepest sites (DCA2 > 3.5), whereas Ericaceae, Asteraceae, and Labiatae towards the less sloping sites (DCA2 < 2); while Rosaceae and Ericacea had the broadest niche widths (2t = 3.4 and 2.7, respectively).

Shrub Functional Group Responses along Coenoclines
Considering the geographical distribution, only Mediterranean shrubs species did not shows a unimodal response curve along DCA1 coenocline of increasing aridity. They showed an increasing trend below maximum attainable response (HOF model III; Figure 3A) towards the south (DCA1 right-end). The shrubs classified as Atlantic, Atlantic-Mediterranean, and endemic of Iberian Peninsula showed skewed unimodal response curves along DCA1 coenocline (model V; Figure 3A; Table 2). Atlantic shrubs with optimum in the northern mountains (DCA1 < 2) had the broadest niche (2t = 4.4), whereas Atlantic-Mediterranean, and endemic of Iberian Peninsula shrubs with narrower niches (2t = 1.7 and 2t = 1.1, respectively) had the optima in the transitional area between the Atlantic and Mediterranean conditions (DCA1 ≈ 2).
When biotype or Raunkiaer s life-form was considered, both deciduous and perennial or phanerophytes and chamerophytes shrubs showed skewed unimodal response curves along DCA1 coenocline (model V; Figure 3B,C; Table 2). Phanerophytes and deciduous shrubs showed close niche widths (2t = 3.0 and 2.6, respectively) with the optima in the northern mountains (DCA1 < 2), whereas perennial and chamaephytes with the same broadest niche width (2t = 4.0) had optima in the transitional area between the Atlantic and Mediterranean conditions (DCA1 ≈ 2).
In the case of shrub dispersal modes, only autochory and barochory did not display unimodal response curves along DCA1 coenocline ( Figure 3D): autochory showed a HOF model II with a decreasing trend towards the south (DCA1 right-end), while barochory showed HOF model III with an increasing trend below maximum attainable response towards the south (DCA1 right-end). Anemochory and zoochory had the optima in the transitional area between the Atlantic and Mediterranean conditions (DCA1 ≈ 2), whereas 'other' had the optimum in the right end of DCA1 coenocline, i.e., towards the south; the three dispersal modes had similar and intermediate niche widths (2t = 2.0-3.3; Table 2).
Both categories of regeneration method displayed skewed unimodal response curves along DCA1 coenocline (model V; Figure 3E; Table 2): resprouters with the optimum in the northern mountains (DCA1 < 2) and germinators with the optimum in the middle part of the gradient (DCA1 ≈ 4); both with narrow niche widths (2t < 1.5).
Along DCA2 coenocline of increasing steepness, and considering shrub geographical distribution, only Atlantic shrubs did not show unimodal response curves but a HOF model III with an increasing trend ( Figure 3F) towards the steepest sites (DCA2 right-end). Shrubs classified as the Mediterranean, Atlantic-Mediterranean, and endemic of the Iberian Peninsula showed skewed unimodal response curves (model V; Table 2) with optima in the steepest sites (DCA2 < 2); the Atlantic-Mediterranean shrubs had the broadest niche width (2t = 3.9).
When biotype or Raunkiaer s life-form were considered ( Figure 3G,H), both deciduous and perennial, as well as chamerophytes showed skewed unimodal response curves along DCA2 coenocline (model V; Table 2) with close niche widths (2t = 1.5-1.8). Chamerophytes and perennial shrubs had the optima in the steepest sites (DCA2 < 2), whereas deciduous shrubs had the optimum in the middle part of the slope gradient (DCA2 = 4.6). Phanerophytes showed a HOF model II with an increasing trend towards the steepest sites (DCA2 right-end).   In the case of shrub dispersal modes ( Figure 3I), anemochory showed a HOF model III with an increasing trend towards the steepest sites (DCA2 right-end). The remainder groups showed unimodal response curves (model V; Table 2) along DCA2 coenocline. Authochory, barochory, and 'other' with very narrow niche widths (2t < 1) had optima in the steepest sites (DCA2 left-end), whereas zoochory with the broadest niche width (2t = 4.6) had the optimum towards the less slopping sites (DCA2 right-end).
Both categories of regeneration method showed skewed unimodal response curves along DCA2 coenocline (model V; Figure 3J; Table 2). Resprouters with the optimum in the steepest sites (DCA2 left-end), and germinators with the optimum towards the less slopping sites (DCA2 right-end); both with narrow niche widths (2t < 1.5). The height of the response (h), i.e., probability of occurrence, was very low along both coenoclines, particularly for taxonomic groups (Table 1). Only Labiatae and Ericaceae along DCA1, and Betulaceae along DCA2 showed probabilities of occurrence above 7%. Probabilities below 20% were predicted for most of the other shrub functional groups along both coenoclines (Table 2).

Discussion
Overall, our findings suggest that understory shrub main families and functional groups respond to the same regional complex environmental gradients as tree species (i.e., altitude temperature, and steepness) since most groups showed unimodal trends to both gradients (coenoclines). However, there were some families ( Figure 2) and functional groups (Figure 3) that showed increasing and decreasing trends along the coenoclines, suggesting that the optimal niches for these groups are out of the length of those gradients. In any case, the description of niche width along gradients for most important shrub families and functional groups is fundamental for two facets: i) to reduce the complexity of the successional process enabling us to understand it better and being able to use sciencebased management for biodiversity and ecosystem services conservation, and ii) to develop sustainable forest practices. Identifying such patterns in the field is crucial to advance ecological forest knowledge taking into consideration not only the tree component but also shrub diversity and function [51].

Distribution Patterns of Main Shrub Families along Coenoclines
The distribution pattern of the main shrub families along the increasing aridity coenocline (DCA1) reflected the recognized environmental gradient from Atlantic temperate areas to the Mediterranean continental climate of the Castilian plateau [44]. The shrub families dominance changes through the first conenocline from the understory shrubby vegetation belonging to Aquifoliaceae, Thymelaceae, and Betulaceae families in the more Atlantic, cold and moist forests of the northern mountains and valleys (DCA1 left-end; Figure 2A) to Ericaceae, Rosaceae, and Cistaceae families in the less moist forests of pre-mountainous range (DCA1 centre), and the more heliophilous shrubby understory with Labiatae and Asteraceae families characteristic of the more xeric Mediterranean forests on limestone and gypsum soils in the south (DCA1 right-end). However, there were two families, Fabaceae and Caprifoliaceae, that showed constant cover values along the main aridity coenocline. It is well known that Fabaceae grow in very different aridity requirements, which favors their ability to inhabit different ecosystems [8]. Similarly, Caprifoliaceae is a family with genera like Sambucus, Viburnum, and Lonicera widely distributed in the northern hemisphere adapted to wide niche requirements [8].
The deciduous broadleaf forests of Atlantic northern areas dominated by Fagus sylvatica and Quercus petraea are associated with an understory of Aquifoliaceae, Ericaceae, and Thymelaceae species such as Daphne laureola L., Erica arborea L., Vaccinium myrtillus L. However, in mountain areas the forests of the moistest valleys are composed of Rosaceae and Betulaceae species mixed with Sorbus aucuparia L. and Ilex aquifolium L. under Betula pubescens Ehrh. As the aridity gradient increased (i.e., moving southwards), more sclerophyllous tree species appeared like Quercus pyrenaica, which has an understory dominated by Ericaceae (e.g., Calluna vulgaris (L.) Hull, Daboecia cantabrica (Huds.) K.Koch, Erica cinerea L., Erica scoparia L.) that constitute the substitution shrubby understory of more transitionalhumid woodlands. In turn, conifer woodlands (Pinus sylvestris, Pinus uncinata Mill. ex Mirb.) in the cold and rainy mountains of the north, shrub-understory is dominated by Ericaceae family with Erica arborea arborea as the main species. Finally, more xerophytic and thermophilic shrub families replaced this understory when moving southwards, such as Labiatae and Asteraceae (e.g., Lavandula spp. and Helichrysum stoechas (L.) Moench) dominating under Quercus ilex subsp. ballota, Quercus faginea, and Juniperus thurifera forests, whereas in the limestone moors of the south, under Pinus halepensis, Pinus pinea, and Cupressuss sempervirens plantations, heliophyllous and summer drought shrub families dominate, such as Fabaceae and Asteraceae (e.g Prunus spp., Dorycnium pentaphyllum Scop., Santolina rosmarinifolia L., and Ononis tridentata L.).
Concerning the slope gradient coenocline (steepness, DCA2), our results showed a turnover of shrub families from those dominating on the limestone moors and flat xeric forest lands where Ericaceare, Labiatae, Asteraceae, and Cistaceae cover the shady soils of Pinus plantations such as Pinus pinaster and Pinus nigra and Quercus ilex stands (DCA2 left-end), to more steeply sloping sites occupied by Betula spp. mixed forest with Aquifoliaceae and Rosaceae on the Atlantic mountain range. Finally, on locations with large rocky outcrops and higher steepness an understory of Oleaceae, Fabaceae, Carpifoliaceae, and Thymelaceae dominates (DCA2 right-end).
Unsurprisingly, the most frequent HOF models for shrub families were unimodal response curves for both gradients (50% of HOF modeled families), as previously found for trees along the same coenoclines (87% of HOF modeled species [41]). Considering the fairly length and complexity of the environmental gradient selected, it seems that there is a large amount of compositional turnover in the vegetation along the north-south environmental gradient and steepness gradients analyzed [41]. These results agree with Økland [52] who stated that the relative frequency of monotonous (HOF models II or III) and unimodal (HOF models IV or V) response curves shift in favor of the latter when the amount of compositional turnover along with gradients increases. In addition, when compositional turnover increased also indeterminate response curves (HOF model I) are expected to be less frequent [52,53]. Our results agree with this expectation since there is a low percentage of species showing monotone response curves (only 20% of HOF modeled families for north-south gradient), as it was also found in Danish forests [49]. It seems that unimodal curves are very more common among plant species along long-gradients with high species turnover [53], thus the shrub response shapes might be mainly determined by gradient properties such as length.
Although there was a greater percentage of symmetric unimodal vs. skewed response curves (80 vs. 20% along DCA1 and DCA2; model IV vs. model V), like expected for long length gradients [53], only a few shrub families showed distinct skew curves (Aquifoliaceae along DCA1 and Betulaceae along DCA2). However, only Aquifoliaceae have its optima near the gradient endpoints according to Rydgren et al. [53], while Betulaceae has its optima in an intermediate position. Additionally, even selecting long gradients a substantial portion of species having truncated realized response curves is probable to identify since their optima will fall outside the sampled part of the gradient [54]. This might be the reason why the response curves shape of Thymelaceae, Oleaceae, or Asteraceae along the north-south gradient (DCA1) and Oleaceae, Thymelaceae, Fabaceae, and Cistaceae along the steepness gradient (DCA2), which showed HOF model V, do not correspond with expected unimodal shape, probably because their niche optima fall outside the environmental gradient associated to this coenocline. This is particularly interesting in the north-south gradient where Atlantic shrub families such as Thymelaceae and Oleaceae have optima to the left endpoint and more marked Mediterranean species such as Asteraceae have it to the right endpoint.

Unimodal Response of Shrub Functional Groups to Environmental Gradients
As expected, the shrub functional groups dominance followed the main tree species trends [41] with a progression along the aridity coenocline from Atlantic high mountain deciduous shrubs with zoochorous dispersion (Sorbus aucuparia and Vaccinium myrtillus), or perennial shrubs (Calluna vulgaris, Daboecia cantabrica, and Erica cinerea) in the north (DCA1 left-end) to more Mediterranean perennial shrubs (Thymus spp., Lavandula spp., Dorycnium pentaphyllum, Santolina rosmarinifolia, and Ononis tridentata), or endemic of the Iberian Peninsula (Chamaespartium tridentatum (L). P.E. Gibbs) towards the south (DCA1 right-end). This distribution pattern might be caused by the different niche requirements of the shrubs described across the coenocline [8]. These results suggest that the shrub community assembly might be ruled by functional differentiation of shrub species [55].

Shrub Optima and Niche Widths across the Main Coenoclines
The lower probability of occurrence of shrub species along both coenoclines, in comparison with tree species (see [41]), might be caused by a lack of exhaustive shrub sampling during the 3SNFI. Here, the main shrub species per field plot were only recorded, i.e., those with a cover above 2%. At the same time, shrub species identification in the 3SNFI is limited to a predefined list of 169 taxons (125 species, 42 genera, and 2 subfamilies), and only individual species are differentiated if they are considered frequent or important enough. It should be noted also as a limitation that shrubs should be successfully identified by the field crews, while the rest of the species are grouped mostly at the genus level [27].
Most shrub species groups with the highest probability of occurrence along both coenoclines (30% > h > 10%) have their optima towards the ends of the ecological gradients, but they do not have the highest niche widths, which is opposite to the outputs when tree species were considered along the same coenoclines [41]. In contrast, most shrub species groups with wide niches along both coenoclines (2t > 1.5) have their optima in the middle part of the gradients.
Rosaceae and Cistaceae or Rosaceae and Ericacea along DCA1 and DCA2, respectively, have particularly wide niches (2t > 2.7). Families occupying environments with sharp contrast and thus present in diverse ecological conditions have broader niches [49]; for example, Rosa spp. in limestone moors and mountains areas, Cytisus scoparius in detrital moors and mountains, and Calluna vulgaris in both forest gaps and dense forests on slightly sloping lands that is a good indicator of soil impoverished conditions.
The narrowest niches found in other families along both coenoclines may be explained by their preference for highly specialized habitats; for example, in mountain areas, Betulaceae and Aquifoliaceae are part of the Atlantic deciduous broadleaf forests characterized by species such as Betula pubescens, Fagus sylvatica, and Quercus petraea. Similarly, in the more temperate and Mediterranean forests of the southern part of the gradient, Thymelaeaceae and Oleaceae are abundant showing the narrowest niches.
Concerning the main taxonomic groups, the Cistaceae family, considered as indicators of Mediterranean conditions [56], has a broad niche width along DCA1, with optimum to the middle-north part of the gradient, marking the division between the Atlantic and Mediterranean conditions. Cistaceae family also shows a preference for less steeped sites (left end of DCA2 coenocline) as well as Labiatae and Asteraceae taxonomic groups. Fabaceae (with constant cover along DCA1) and Rosaceae (with broadest niche widths along both coenoclines) are more generalists [57] and distributed through almost the whole DCA1 and DCA2 coenoclines, although their cover increase towards the steepest sites, in both Atlantic and Mediterranean locations.
As expected, niche locations of certain functional groups of shrubs are close in DCA1 (zoochory and anemochory or geographic range, Table 1; [55]). Similarly, there is a considerable overlap of ranges between shrubs geographical distribution along the first coenocline, indicating that the geographical distribution ranges do not show edge limits from north to south of the gradient. The rather balanced relations between the forest shrub species, indicating species co-existence [49], are common, e.g., Thymelaeaceae and Fabaceae in limestone moors, Ericaceae in detrital moors or the mountains. These results seem to support the idea that shrub families (Cistaceae, Rosaceae, Ericaceae) or functional groups (e.g., phanerophytes and deciduous shrubs) with niches located in the central DCA1 coenocline can move northwards (DCA1 left end) in response to predicted climate change scenarios [30]. In contrast, shrubs growing at the mountainous edge of the gradient (e.g., Aquifoliaceae, Betulaceae) have more specialized niches and can only move in altitude [30].
When the coenocline related with steepness was considered, we observe that niche location of several functional groups (e.g., perennial vs. deciduous, chamaephytes vs. phanerophytes, resprouters vs. germinators) were very distant with optima around lower values for perennial, chamaephytes, and resprouters shrubs, whereas deciduous, phanerophytes, and germinator shrubs are mainly located at higher optima values, indicating that shrubs conforming these groups are specialist of specific habitat positions (steepness). In contrast, Atlantic shrubs, with a wide distribution along the steepness gradient, and being more generalists showed a strong asymmetry in the overlap indicating some preferences towards higher steeped sites [8].

Conclusions
The topographic-climatic gradient (north-south primary coenocline) and the slope gradient (secondary coenocline) are strong environmental gradients that shape, not only the tree responses but also the shrub families and functional groups responses. Therefore, along these two main coenoclines identified it is possible to relate dominant tree species with their understory shrubs. Similarly, our results also provide evidence for the generality of unimodal shrub family response curves, when the gradient selected includes a sufficient amount of compositional turnover. In the same way, shrub groups inhabiting environments with sharp contrast or transitional environments (i.e., occupying different morpho-structural units) have the broadest niches, similarly to those described for trees. In contrast, those species occupying localized habitats showed the narrowest niches. These results have significant implications for forest management since shrub families and groups with the narrowest niches are the most sensitive, with a potential greater loss of climate niche space considering the present climate change projections. At the same time, the species with broadest niches can be used as pioneer shrubs to increase heterogeneity in micro-environmental conditions enhancing forest expansion. Finally, the combination and analysis of forest inventory data that have been collected at different dates will provide tools to practitioners for the identification of changes in distribution patterns and niche sizes for the species of interest, being possible to relate, easily, species distribution changes with environmental or management changes. Using this family or functional approach for the shrub layer, we have solved the main limitation in the 3SNFI regarding taxonomic identification of shrub species, being a simple methodology to give added value for forest planning to the shrub information collected in the forest inventories.    [58]. (b) Biotype and Raunkiaer's life-forms classification following [59]. (c) Dispersal mode according to LEDA traitbase [60] and D-3 database [61]. (d) Regeneration method according to [62][63][64][65].