Pastoral and woodcutting activities drive Cedrus atlantica Mediterranean forest structure in the Moroccan Middle Atlas

. Human activities are historical ecological drivers, and we need to better understand their effects on ecosystems. In particular, they have been very important in the shaping of the Mediterranean biodiversity hotspot. Researchers and managers nonethe- less lack knowledge concerning the impacts of their combinations and their current intensity on the structure of forest ecosystems of the southern part of the Mediterranean basin. In this study, we have develped a new methodology in order to understand the impacts of combined pastoral and woodcutting activities on the forest structure of the still ill-described but ecologically and economically important Moroccan Middle Atlas cedar forests. In a 40 000 ha forest, we chose 103 sites and sampled human activities through proxies and forest structures through circumference and vertical structures. A typology of sites yielded four human activity types: dominant pastoral activities, dominant oak cutting or cedar cutting activities, and an intermediate mid-disturbance type. This typology did not depend on altitude or substrate, confirming that the ecosystem structures linked to the different types depend more on human activities than on main environmental parameters. Pastoral activities modified forests the most, converting them to parklands with reduced canopies and low dynamics but high tree maturation. Woodcutting activities induced gap dynamics, favoring Cedrus atlantica in favorable environmental conditions and Quercus ilex otherwise, while they affected vertical structure depending on the local environment and competition for light and soil resources. Moderately disturbed stands showed forest maturation with low competition for light. Unlike previous studies, we found no evidence of a general degradation of cedar forests due to local human activities. However, cedar logging has reduced standing basal area regionally and one third of the sites may have vulnerable cedar populations due to pastoral activities and to unfavorable environmental conditions. These results can direct future research and management needs for a better protection of Mediterranean forests and parklands and their biodiversity, although to be effective such efforts must also partner with sociogeographical studies.


IntroductIon
Scientists are now realizing how widespread human activities have been in ecosystems and the need to take them into account to understand their functioning (Willis et al. 2004), especially in the Mediterranean basin. In this anciently occupied region, humans have tapped the ecosystems for resources for millennia, modifying ecosystem characteristics (Grove andRackham 2001, Quézel andMédail 2003) while maintaining high levels of biodiversity (Blondel 2006). However, this biodiversity is threatened by the speed of both climate and anthropogenic changes (Quézel and Médail 2003). On the northern shore, biodiversity is menaced by woody encroachment and habitat fragmentation (Scarascia-Mugnozza et al. 2000). Conversely, on the southern shore, demand for natural resources is ever growing and human activities inside forests are intensifying (Scarascia-Mugnozza et al. 2000). In this context, both researchers and managers critically lack knowledge about the impacts of human disturbances on Mediterranean forests (Scarascia-Mugnozza et al. 2000) and especially their combined effects (Wisdom et al. 2006). This lack of knowledge is especially true on the southern shore, where empirical data is scarce. Understanding these forests' functioning could (1) give insights considering pressure levels necessary to maintain biodiversity in the Abstract.
Human activities are historical ecological drivers, and we need to better understand their effects on ecosystems. In particular, they have been very important in the shaping of the Mediterranean biodiversity hotspot. Researchers and managers nonetheless lack knowledge concerning the impacts of their combinations and their current intensity on the structure of forest ecosystems of the southern part of the Mediterranean basin. In this study, we have develped a new methodology in order to understand the impacts of combined pastoral and woodcutting activities on the forest structure of the still ill-described but ecologically and economically important Moroccan Middle Atlas cedar forests. In a 40 000 ha forest, we chose 103 sites and sampled human activities through proxies and forest structures through circumference and vertical structures. A typology of sites yielded four human activity types: dominant pastoral activities, dominant oak cutting or cedar cutting activities, and an intermediate mid-disturbance type. This typology did not depend on altitude or substrate, confirming that the ecosystem structures linked to the different types depend more on human activities than on main environmental parameters. Pastoral activities modified forests the most, converting them to parklands with reduced canopies and low dynamics but high tree maturation. Woodcutting activities induced gap dynamics, favoring Cedrus atlantica in favorable environmental conditions and Quercus ilex otherwise, while they affected vertical structure depending on the local environment and competition for light and soil resources. Moderately disturbed stands showed forest maturation with low competition for light. Unlike previous studies, we found no evidence of a general degradation of cedar forests due to local human activities. However, cedar logging has reduced standing basal area regionally and one third of the sites may have vulnerable cedar populations due to pastoral activities and to unfavorable environmental conditions. These results can direct future research and management needs for a better protection of Mediterranean forests and parklands and their biodiversity, although to be effective such efforts must also partner with sociogeographical studies. northern part; (2) help conservation both locally and in vicariant ecosystems.
A hotspot within the Mediterranean biodiversity hotspot (Mittermeier et al. 1998, Médail andDiadema 2006), the Moroccan Middle Atlas cedar forests are ecologically important and play a major economical role regionally and nationally. Ecologically, they are North Africa's largest continuous forest ecosystems and host many endemics including three emblematic species: the Atlas cedar (Cedrus atlantica), the thuriferous juniper (Juniperus thurifera), and the Barbary macaque (Maccacus sylvestris; Quézel and Médail 2003). Economically, they are at the heart of pastoral and woodcutting activities and, as such, the base of the ovine and craft-wood value chains. Many witnesses consider these forests to be threatened (BCEOM-SECA 1996, Quézel and Médail 2003, M'Hirit 2006, but such claims are contested (Davis 2005). Up to now, some studies have worked locally on the impacts of woodcutting activities (Linares et al. 2011a, Navarro-Cerrillo et al. 2013. But pastoral activities (grazing and pruning of trees) were never considered in studies on their dynamics despite their critical role for tree recruitment in Mediterranean areas (Plieninger et al. 2003). A research problem in many southern shore areas is that no controlled conditions can be obtained naturally or experimentally and reliable data about activities are scarce. Researchers therefore need to devise new methodologies that can effectively characterize human activities locally and determine their effects on ecosystems.
In this study, we sought to understand the combined impacts of pastoral and woodcutting activities on cedar forests, developing a new methodology and new empirical and theoretical insights that can be used to understand and manage forests locally and elsewhere in the Mediterranean basin. To achieve this, we sampled the biomarkers (stumps, pruning scars, herbivory induced deformation on shrubs) human activities leave in the forests (Étienne andRigolot 2001, Furukawa et al. 2011). In the same places, we sampled forest structures through size class distributions, giving insights into dynamics taking place in ecosystems (Pulido et al. 2001), and vertical structures, yet poorly documented in Mediterranean forests even though they are very important to the understory conditions regulating recruitment and biodiversity (Van Pelt and Franklin 2000). We hypothesized that (1) the variability of human activities, their intensities, and likely their combination could be clustered into distinct activity types. Based on previous studies (Pujos 1964), we assumed that natural physical environmental parameters, mainly altitude, have a strong effect on forest structure and stock. Therefore, we also hypothesized that (2) the variability of human activities is linked to altitude, this variable summing up the main environmental conditions; and (3) the activity types are linked to different forest structures, which are linked to different dynamics. Considering some alarming data about cedar forests (BCEOM-SECA 1996, Quézel and Médail 2003, M'Hirit 2006, we finally hypothesized that (4) some human activity types induce local cedar population vulnerability, assessed through a weak (rather flat instead of J-shaped) dynamics.

Study area and sampling sites
The Moroccan Middle Atlas hosts the largest C. atlantica forests, with a total of nearly 80 000 ha (Appendix S1). This species grows in altitudes between 1500 and 2300 m (Quézel and Médail 2003) in a humid to subhumid mediterranean climate with high inter-annual variability (mean January minimum temperature in Ifrane at 1600 m, 2.7 • C; mean July maximum temperature, 28.7 • C; mean annual precipitation between 1500 and 500 mm decreasing from West to east; Nègre 1952, M'Hirit 2006. Substrates are mainly carbonated, but cedar forests also grow on acidic substrates (Quézel and Médail 2003). C. atlantica grows in mixed stands with Quercus ilex sp. rotundifolia, and several minor companion species (Quézel and Médail 2003).
Considering how widespread human activities are in the Middle Atlas and the difficulty to obtain a priori information on their local intensities, we chose to measure human activities and ecosystem structure response in situ. In the spring of 2013, 103 sites of 400 m 2 were sampled on carbonated substrates in four administrative zones covering approximately 40 000 ha of cedar forests. Patterns of human activities and forest structures were sampled semi-systematically along transects constituting a random grid covering the entire study region (Appendix S1). Transects ran west to east along the regional increase in altitude. For each administrative zone, nine replicate sites were sampled for each of three altitude classes. One zone had only the highest altitude class, so 10 sites were sampled for that class. In the northern zone, 13 sites sampled in a preliminary study were also used.

Ecosystem structure sampling
For each site, tree size class structure was obtained through circumference measurement at 1.3 m of all stems of circumference larger than 25 cm (3394 trees sampled, 26% C. atlantica, 72% Q. ilex, and 2% other species). This data yielded stand parameters: specific or total basal area and stem numbers as well as ratios. Resembling Plieninger et al. (2011), current C. atlantica recruitment was measured: saplings, stems that had survived their first year and with circumference at 1.3 m under 25 cm, were counted. C. atlantica seedlings and seeds were counted in five 2 × 2 m plots evenly spaced inside the site. Seeds were counted using an ordinal variable (four categories: absent, 1-10 seeds, 10-30 seeds, 30 seeds and more).
For each tree, the crown was reconstructed using the maximum crown height measured by a Haglöf Vertex (Långsele, Sweden) and minimum crown height and radius estimations (all done by a single observer so the bias, if it existed, was similar throughout the sampling). Vertical structure was reconstructed based on the "crown shell method" described in Parker et al. (2004) as it was considered most efficient to link prunings directly to the affected individual. For each tree, the crown shell method models a crown shape by fitting a predetermined geometrical shape to crown measurements (Parker et al. 2004). Unlike other studies with other species in Morocco (Montès 1999), our study included too many trees for a precise accountability of individual shapes. As trees mostly had elongated shapes and no precise predetermined shape could be applied, cylindrical shells were used. Vertical stand structure was then recomposed by summing the individual crowns. We thus derived vertical structure parameters: heights, canopy length and volume, and canopy area projection (Table 1).

Human activity sampling and disturbance indexes
Few data are available regarding the main human forest perturbations: pastoral and woodcutting activities. Pastoral activities are of two types: forest understory grazing by herds of 10-500 sheep and goats, and prunings for fodder when grass is sparse to avoid buying feed (BCEOM-SECA 1996). Officially, grazing is limited to familial livestock (informally considered to be 50-60 heads), and pruning is illegal. This results in the near absence of official or statistical data regarding these activities. In the region, woodcutting activities are led for different purposes (firewood and charcoal, construction poles, high value craftswood). Again, as there is much illegal cutting, Forest Service data cannot be used as a surrogate sampling basis.
The biomarkers left by human activities were therefore sampled and converted to seven numerical indexes table 1. Table of human disturbance indexes and their response stand and vertical structure parameters (mean ± SD).
interpretable as perturbations (noted 1-7). Grazing intensity was estimated using the deformation livestock induce on live Q. ilex shrubs (Étienne and Rigolot 2001; five categories: unbrowsed, browsed but undeformed, browsed and deformed, severe deformation, leaves still present, and severe deformation with disappearance of most leaves), and was converted into an ordinal grazing index (1). Pruning intensity was measured through scar counts on each tree. As cutting a few branches does not have the same effect on a small tree than on a large one, pruning was converted into an intensity. Pruning disturbance indexes were the mean scar counts standardized by canopy volume of each tree for C. atlantica and for Q. ilex (2 and 3, respectively). Woodcutting activities were sampled by species and circumference of every stump larger than 10 cm (Furukawa et al. 2011). The woodcutting induced disturbance depends on the removed to remaining biomass ratio. Stump circumferences were therefore converted into basal area of stumps and standardized per species by standing basal area (4 and 5 for C. atlantica and Q. ilex, respectively). The type of disturbance is not the same if small trees or big trees are removed for the same proportion of removed basal area. Numbers of C. atlantica and Q. ilex stumps were therefore the last two woodcutting indexes used (6 and 7 for C. atlantica and Q. ilex, respectively). These were correlated to stump basal areas, also giving an idea of the absolute cutting that took place. The measured stumps could be dated neither through the unreliable forest administration data, nor through dendrochronology methods due to their high numbers. A subsampling analysis is ongoing (B. Brossier, personal communication), and preliminary results show stumps record disturbances that are relatively recent: over 50 yr for C. atlantica, less but difficult to determine for Q. ilex.
Apart from pastoral and woodcutting activities, humans extract non-timber forest products in the forest understory. We did not sample them as they are led irregularly throughout the region and they have low impact on forest structure. Currently in Morrocan cedar forests, human-induced fire was minor and we removed the only site that had been burned from analyses.

Statistics
All statistics were computed with R (R Core Team 2014).
To test hypothesis 1, a stand typology according to human activities was obtained by running a cluster analysis with the partitioning around medoids method on the axes of a principal components analysis (PCA) computed from all human disturbance indexes at all sites. The cluster number was chosen to maximize silhouette profiles (Kaufman and Rousseeuw 2009). The obtained types were compared through several Kruskall-Wallis tests on human disturbance indexes. To test hypothesis 2 and the effect of altitude on the partitioning, Mann-Kendal trend tests (Hamed and Rao 1998) were run between site altitude and position on PCA axes. The potential effect of other environmental parameters was tested in Appendix S2.
To test hypothesis 3, and consequently hypothesis 4, we used stand parameters, size class structures of live trees and stumps combined with recruitment parameters and vertical structures. Kruskall-Wallis tests were run on stand and vertical structure parameters, using C. atlantica and Q. ilex only because other species were not abundant enough for consistent analyses, even as indicator species. Significant differences were also tested among types for C. atlantica regenerations (saplings, seedlings, and seeds), while correlations between regeneration stages were tested using the Pearson correlation test (except for the ordinal seed numbers for which the Spearman correlation test was applied).
Size class structures were obtained for live trees and stumps and compared among human activity types and to the regional population, obtained by averaging parameters of all sites. Circumferences were binned into 15-cm classes to obtain tree and stump size class structures. To assess the difference between each type size class structure and the averaged one for the region, relative anomalies were constructed (Population value−Type value/ Population value) and curves smoothed using the loess algorithm. Positions of anomaly curves to the zero show whether the type and the regional population have different numbers of stems or of stumps (absolute woodcutting). With positive anomalies, more stems or stumps are present in the type than in the population. Within a type, the position between stump and stem anomaly curves show whether type and regional populations have different stump to stem ratios (proportional woodcutting). Indeed, if the stump curve is above the live tree curve (high proportional woodcutting), more trees have been cut per available resource in the type than in the regional population. In Appendix S3, stump size class structures were also obtained for comparison.
Observed vertical structures were compared among human activity types using Kruskall-Wallis tests for each 0.5 m height interval. To further test the impacts of human activities on canopies, observed structures of each type and of the regional population (total number of pooled trees) were compared to reconstructed unpruned vertical structures. To proceed, for each type, we used the same size class structure and randomly selected unpruned trees from the regional population (the 736 unpruned C. atlantica and 1870 unpruned Q. ilex measured in this study). For statistical comparison, 1000 reconstructed vertical structures were obtained through the following random resampling method: each pruned or unpruned tree of each activity type was randomly reassigned the crown parameters of a tree in the same size class from the regional dataset comprising only unpruned trees. If no unpruned tree existed in the size class, we assigned the mean class canopy size obtained by a cubic regression on all unpruned trees. For each 0.5 m height interval, if the observed type structure was not within the 95% confidence interval shell of the 1000 reconstructed vertical structures, results were considered to be significant. This innovative method enabled us to test changes due to prunings in heavily pruned areas and to test whether for a same size class structure the human activity type's vertical structure differed from the vertical structure expected with random unpruned trees of the regional population. Results could give an over-representation of unusual crown parameters in the larger size classes where unpruned trees are few. However, these trees were rare and the effect was assumed to be weak. Moreover, considering the difficulty of controlling field conditions in our region, this approach was the best available to reconstruct stands unaffected by human activity types.

Human disturbance typology
Based on silhouette profiles, the cluster analysis provided six human activity types. Two were too small for subsequent statistical analyses and were merged into the closest ones resulting from the PCA (Fig. 1). Confirming hypothesis 1, the differences among the four remaining human disturbance types stood out using the Kruskall-Wallis tests (Table 1A). The first, including 17 sites, was the pastoralism type (PA) as it had significantly larger values of grazing intensity and globally larger pruning intensities on C. atlantica and Q. ilex. The second, with 12 sites, was the oak cutting type (OC). It had significantly larger numbers of Q. ilex stumps and larger Q. ilex cutting intensity than all other types. The third, with 15 sites, was the cedar cutting type (CC). It had significantly larger numbers of C. atlantica stumps than all other types and a significantly larger cedar cutting intensity than PA and than the last type. The last type, with 59 sites, was the mid-disturbance type (MD). It had significantly smaller values of pastoral indexes than PA and smaller cutting indexes than CC and OC but larger grazing intensity than CC. Despite significant differences in human disturbances, few significant values were found in usual forest stand parameters (Table 1B). OC differed from MD with smaller total basal area and C. atlantica basal area, and PA stands differed from CC by a smaller number of C. atlantica. This relative absence of significant differences confirmed the importance of looking further into size class and canopy structures.
FIg. 1. Typology of study sites according to human disturbances projected on first PCA plane (47% of inertia). The four ellipses represent the four human activity types with associated forest sites (dots) and prevailing human disturbance indexes (arrows). Abbreviations are RelCedarStArea, relative Cedrus atlantica stump area index; RelOakStArea, relative Quercus ilex stump area index; NOakStumps, number of Q. ilex stumps; NCedarStumps, number of C. atlantica stumps.

Effect of altitude on human activities
Disconfirming hypothesis 2, there were no significant trends between altitude and PCA axes. Similarly, the projection of sites on a PCA of other environmental variables gave no clear differences among types (Appendix S2). Only PA type and CC type differed significantly, PA type being in hotter and less continental places.

Structures and dynamics linked to human activity types
Size class distributions At the regional scale, both C. atlantica and Q. ilex size class distributions were inverse J-shaped (Fig. 2). This pattern was true for both species for all types with one exception: the C. atlantica distribution of PA type was flat, showing low dynamics. OC type C. atlantica stands also had reduced dynamics. These two types may therefore have vulnerable C. atlantica populations confirming hypothesis 4.
Anomalies for both species of the PA type showed similar trends with low live tree abundance under 400 cm circumference for C. atlantica and 200 cm for Q. ilex, and high abundance at large sizes, and with low absolute and proportional woodcutting (Fig. 2). Anomalies of the OC type showed a low abundance and low C. atlantica dynamic with a high proportional woodcutting, whereas the Q. ilex population was dynamic with abundant trees smaller than 100 cm and low abundance at larger sizes. Absolute and proportional woodcutting were high. Anomalies of both species of the MD type showed the same trends with a high abundance of large trees (from 100 cm to 500 cm circumference for C. atlantica, and Q. ilex larger than 80 cm circumference), while absolute and proportional woodcutting were low (Fig. 2). Anomalies of the CC type showed inverse specific profiles with a very dynamic C. atlantica population with a low abundance of trees larger than 250 cm in circumference and high woodcutting, whereas the Q. ilex population dwindled to zero for circumferences larger than 150 cm and had low cutting.
Apprehension of current dynamics was given by differences in present recruitment (including saplings, seedlings, and seeds). No significant differences existed for recruitment parameters among the four types (Table 2). However, the number of saplings increased from PA to CC. Moreover, there were positive correlations between sapling numbers and categorical seed numbers (Spearman's P < 0.001, ρ = 0.32) and between seedlings numbers and categorical seed numbers (Spearman's P < 0.001, ρ = 0.57), but there was no correlation between sapling numbers and seedling numbers. The proportion of sites with seeds was always higher than the proportion of sites with saplings, which was higher than the proportion of seedlings.

Vertical structures
Observed vertical structures showed that species stratification patterns were similar among types (Fig. 3, left panels). Logically, the canopy was dominated in height and volume by C. atlantica trees with heights up to 45 m. Q. ilex were smaller with heights up to 22 m. They constituted, with small C. atlantica trees, the lower half of the canopy. However, some shape differences exist among types. In MD and CC types, canopy shapes were similar to that of the regional population with wide C. atlantica profiles, widest between 5 and 10 m (Fig. 3, left panels). PA and OC types had narrower C. atlantica profiles. OC type nevertheless maintained a maximum area that was close to the regional population shape with a wide occupancy in lower strata. Profiles were all unimodal, except for the CC type in which C. atlantica had a slightly bimodal shape. The MD type canopy had a significantly larger canopy volume, length and median height than OC type (Table 1C). The only differences found among human disturbance types were for MD observed C. atlantica occupation that was significantly larger than that of OC from 6.5 to 30 m, and that of PA from 7.5 to 28 m (Fig. 3, left panels).
Species patterns between paired observed and theoretical unpruned structures were variable (Fig. 3, paired left and right panels). The PA type showed the most affected profile for both species with a depleted foliage occupation at most heights above 5 m for C. atlantica and above 2.5 m for Q. ilex. Conversely, C. atlantica foliage showed a significant excess at 2.5 m. The OC type was the second most affected with depleted cedar foliage at intermediate heights. CC type had a small depletion around 10 m, corresponding to the bimodal curve's inflection height. Some stands had excentered canopies. OC type oak canopy was higher than expected from the random population, with foliage excess in upper strata and depletion in lower strata. Inversely, MD type cedar canopy and both MD and CC type oak canopies were lower than expected. These differences in size class and vertical structures confirmed hypothesis 3, as activity types could be linked to different forest structures and dynamics.

Four types of forest uses resulted from the two main economic human activities
We found an interaction gradient between pastoral and woodcutting activities. Confirming hypothesis 1, the data was partitioned into four main forest use types (Fig. 1). The type having the most impact on ecosystems was the PA type (representing 17% of sites) where grazing and pruning converted forests to parklands by limiting C. atlantica dynamics and the canopy volume, while favoring large trees. The second in terms of impact was the OC type (12% of sites) where woodcutting led to an undynamic C. atlantica population that seemed to be at the border of its environmental niche and was giving way to a developing Q. ilex population. In the MD type (57% of stands), moderate activities left a maturing forest with low competition for light. Finally, in the CC type (15% of sites), cedar timber cuts and thinnings led to a dynamic C. atlantica cohort previously suppressed by larger trees. It had inverse C. atlantica dynamics than PA with more stems but a similar basal area.
FIg. 2. Size class distributions of C. atlantica and Q. ilex for the regional population and the four human activity types (mean number of stems/ha ± SE). In histograms, gray is trees and white is saplings; solid and dashed lines show type anomalies to regional tree and stump populations, respectively. Positions of the anomaly curves to the zero show whether the type and regional population have different numbers of stems or of stumps (absolute woodcutting). For instance, if the value is −1, type or population abundance is 0, if value is 1 it is twice the population's. The relative positions of the anomaly curves (dashed line above solid line for high proportional woodcutting) show whether more woodcutting took place proportionally to available stock in the type than in the regional population (with an exception for cedar cutting oaks above 150 cm, an artifact due to the smoothing process).
We can state that the structure differences reported among types depended more on human activities than on the expected main environmental parameters. Indeed, rebutting hypothesis 2, altitude was an important factor neither in the PCA structure nor in the resulting types, while the sampling design was only applied on carbonated substrates (the most prevalent in cedar forests). Similarly, there were little differences in environmental parameters among types (Appendix S2). This was a surprise as altitude is considered to be a determinant of forest types (Pujos 1964), which in turn, as a resource, should drive human activities. However, stand parameters did not  differ significantly among activity types, indicating that resource extraction did not directly depend on local forest stock. This demonstration of the impact of human activities on ecosystem structures was made possible thanks to our methodology sampling human activities directly in ecosystems. As expected, this allowed to understand precise patterns of the heterogeneous use of forests that other studies could not completely take into account (Linares et al. 2011b). It also prevented from assuming forest dynamics without an in situ link to human activities as some recent studies have done using forestry data or assessing activity types without verification (Navarro-Cerrillo et al. 2013).

New ecological insights from cedar forests
Regional forest dynamics At the regional population scale, the J-shaped C. atlantica and Q. ilex size class structures could be interpreted as dynamic and healthy populations. This should relativize the catastrophic discourse about cedar forests that many observers tend to have (BCEOM-SECA 1996, Quézel and Médail 2003, M'Hirit 2006, especially considering their dynamics. As C. atlantica is a long-lived species, temporary recruitment in patches is expected to be sufficient to maintain a dynamic population over a 40 000 ha forest (Pujos 1964). However, the regional cedar stump to standing basal area ratio of 0.8 shows woodcutting has had important effects on forest stock. The limiting factor of recruitment was not clear. Parameters suggested a control between the seedling and sapling stages, which were uncorrelated (Table 2). Low sapling abundances were linked to high grazing intensities or low canopy covers, both known to critically affect this stage transition (Lepoutre 1964, Pujos 1964, Quézel and Médail 2003. However, our spring snapshot sampling could not distinguish between the two effects, which are likely interactive.
Confirming hypothesis 3, changing scale and comparing forest stands among human activity types showed strong differences in the associated structures and dynamics.

Parkland ecosystems created by pastoral activities
Cedar forest structure was most affected by pastoral activities, acting like a stress to limit C. atlantica small class survival (Fig. 2) and depleting canopies (Fig. 3). Tree size class distributions indicated PA stands were actually parkland ecosystems dominated by two tree populations with various sensitivities to human activities (Lykke 1998). This pointed to different dynamics at small and large stages.
Considering the flat population dynamics, low but continuous C. atlantica long-term recruitment limited population levels below those of self-thinnings (Coomes et al. 2003, Clark et al. 2012). These population dynamics were likely limited by two stress factors: direct mortality through grazing of seedlings and saplings, and reduced survival to droughts in lower altitudes due to sparse overstory cover (Lepoutre 1964). The latter was confirmed by the location of PA stands in hotter and more continental places than CC type (Appendix S2). Conversely, the J-shaped Q. ilex population indicated self-thinning after low sensitivity of recruitment to pastoral activities. Other studies in Mediterranean parklands reported different Q. ilex spp. grazing related recruitment limitations (Pulido et al. 2001, Plieninger et al. 2011. Their results could be linked to dryer conditions and/or to more intense grazing, while our data could not rule out a recent recruitment limitation that could in the future lead to similar distributions (Plieninger et al. 2011).
On the other hand, pastoral activities were especially compatible with stand maturation. High abundances of large trees of both populations showed low disturbancedriven mortality (Coomes et al. 2003), with low woodcutting in already sparsely stocked areas. However, in the global warming context, this low sensitivity of larger trees to pastoral activities may hide future mortalities (Linares et al. 2011b). Recent decades also had intense prunings and trees might die with a delay.
We showed that pastoral activities affect vertical structure through direct and indirect effects. Expectedly, trees in PA stands had smaller crowns than trees in other types because they were heavily affected by prunings (Fig. 3). For C. atlantica, the effect remained significant well above the species' maximum pruning height (20-25 m), possibly resulting from a bias in our vertical structure reconstruction with cylinders, or hiding a pruning-induced crown growth reduction. Although Q. ilex rebranches (Pulido et al. 2001), this did not compensate for the removed biomass. Finally, the excess in C. atlantica biomass at 2.5 m was unexpected. This might result from a better light availability allowing small trees or rebranching of lower parts of pruned trees to develop better.

Dichotomic dynamics induced by woodcutting activities
OC and CC type size class distributions showed that woodcutting induced gap dynamics whose main recruiting species depended on environmental conditions. Indeed, both types have abundant stumps of all sizes, and CC and OC types have strong J shapes of respectively C. atlantica and Q. ilex (Fig. 2). The most likely scenario is therefore gap dynamics (McCarthy 2001): first, large trees were cut. A release of competition then took place with recruitment and growth of young trees, which were finally thinned. The low grazing intensity in CC compared to MD might have facilitated recruitment, but it was likely not the primary initiator. Indeed, dynamics in the CC type seemed to be a pulse resulting from a temporary lapse in competition: removing the site with the most saplings led the sapling mean to be inferior to the first tree size class.
Although the two types showed cutting of both species, their inversely dynamic Q. ilex and C. atlantica populations (Fig. 2) indicated the resulting gap dynamics depended on micro-site conditions (climate, soils;Nègre 1952, Pujos 1964. Indeed, there were no significant differences on the environmental parameters among types, but OC type tended to be more restricted to low elevation zones with high precipitation while CC type stands were spread out in relatively cool zones (Appendix S2). Although it could be expected that cedar population in OC type would be favored by the cutting of oaks, it was weak both in terms of basal area and canopy volume. Reciprocally, Q. ilex in CC type had a dwindling population. This suggests that these populations were on their respective niches' edges, inducing low C. atlantica dynamics in the OC type that were probably even lower due to depleted canopies in low elevations (Lepoutre 1964) and ongoing global warming (Linares et al. 2011b).
Two recent local-scale studies concluded inverse successional paths after woodcutting (Linares et al. 2011b, Navarro-Cerrillo et al. 2013. Our results showed that at the regional scale, dynamics of both species could take place, although in different sites (OC and CC types). They could however also mature together (MD type). Linares et al.'s (2011b) results would correspond to our OC type. As the latter's sampling did not distinguish basalt and carbonated substrates (Navarro-Cerrillo et al. 2013; pictures in supporting information), their results cannot be interpreted.
Vertical structures of OC and CC types were modified by competition and environmental conditions. Contrary to expectations for thinnings (Aussenac 2000), OC type Q. ilex had higher canopies than average (Fig. 3). The large root systems benefitting the remaining stems (Quézel and Médail 2003) could lead to an exacerbated competition for light and displace crowns upwards (Lines et al. 2012). In contrast, in the C. atlantica population of CC type, the bimodal curve (Fig. 3) indicated that large trees competitively suppress the small ones (McCarthy 2001). Unlike expectations for isolated trees situated in disturbed stands (Giuggiola et al. 2013), the reduced canopy of OC type C. atlantica (Fig. 3) likely resulted from unfavorable environmental conditions. CC type Q. ilex showed a pattern resembling the maturing populations of MD type (Fig. 3).

Low competition for light with low canopies and maturing stands in moderately impacted stands
Forest maturation was important in the MD type (Fig.  2), suggesting a low natural disturbance-driven mortality (Coomes et al. 2003) and showing that woodcutting was likely the main large tree mortality factor at the scale of the forest. Dynamics appeared to be suppressed by competition with larger trees (Pujos 1964, McCarthy 2001. In these stands, low human activities allowed few individuals from some companion species, such as Juniperus thurifera or Quercus canariensis, to attain remarkable circumference and crown sizes.
The vertical profiles of low cut stands (MD type cedar and MD and CC type oak) were likely shaped by the interaction of physiological and ecological constraints due to the Mediterranean climate. Unlike expected high competition for light in low disturbance stands (Lines et al. 2012, Hawthorne et al. 2013, these populations' canopies were displaced downwards (Fig. 3), without link to particular climatic conditions (Appendix S2). This inverted structure could result from relatively open stands with high insolation where larger trees could be less affected by competition for light than by the summer hydric constraint (Ryan et al. 2006). As competition is known to limit water stocks (Aussenac 2000) and increase drought related mortality (Giuggiola et al. 2013, Ruiz-Benito et al. 2013, it could be a further driver of the maturing forests' vertical structure.

Lessons for the Mediterranean basin
Up to now, cedar forests have mainly been studied in terms of "natural" dynamics (Nègre 1952, Lepoutre 1964, Pujos 1964, Quézel and Médail 2003. Our study is a benchmark that could help update this data and direct research and management needs for Mediterranean ecosystems both on the southern and northern shores, especially in systems where one species is preferred for economical/ecological/symbolic reasons and to attain sustainable development goals that must respect both socio-economic priorities and habitat and biodiversity conservation. In the middle Atlas cedar forests, regional forest vulnerability seemed to be limited, with 70% of stands in apparently good health (MD and CC types). However, 30% of sites had low C. atlantica dynamics. Confirming hypothesis 4, these might possibly be vulnerable populations, with 20% being related to intense pastoral activities and 10% to environmental conditions (OC type).
Mediterranean forest dynamics could be vulnerable because of pastoral activities. As these stands corresponded to various environmental conditions, a larger part of the forests could be converted to parklands due to recruitment limitation linked to increased pastoral activities. However, these activities should be maintained not only for socio-economic reasons (BCEOM-SECA 1996), but also because Mediterranean parklands, especially when older (Bergmeier et al. 2010), have a high conservation value for their biodiversity (Grove and Rackham 2001). Pastoral rotations could be prescribed to avoid threatening them by sustained intense pastoralism (Plieninger et al. 2003), but local application has not and will not be easy due to complex social dynamics (Aubert 2010).
We showed that while woodcutting activities can increase the species of interest's dynamics, they can also make populations vulnerable through replacement by another species in the niche edge areas. This highlights that Mediterranean forests are complex with no homogeneous post-disturbance growth. Currently, local forest management aims at favoring gap dynamics and converting MD type forests to CC type to improve their production and reduce senescence (TTOBA 2002). This is sustainable for the persistence of cedar forests in stands that are not on the species' niche edge, but rejuvenating stands might lead to a loss of some ancient Mediterranean forests and their biodiversity (Vallauri et al. 2002, Quézel andMédail 2003).
Several slight sampling limits could prevent generalizing our results. By including many forest stands in a large area, we could not have the precision needed to date all woodcutting or prunings. Data concerning small stems and C. atlantica and Q. ilex recruitment was incomplete and recruitment takes place episodically (Pujos 1964). The limitation stage of recruitment and its reason therefore remain elusive. Some human-induced changes affecting forest structure were not taken into account: past climatic changes have affected cedar forest migration (Cheddadi et al. 2009), and the lower distribution limit is currently moving up with climate change (Rhanem 2011). Even though they were marginal in our region, fires are also known to affect C. atlantica dynamics (Pujos 1964), especially in Algeria. Our study was also incomplete in the apprehension of environmental parameters. Indeed, we only sampled carbonated substrates, while forests on acidic substrates locally and in the Rif mountains are more dynamic (Pujos 1964). Similarly, our climatic range was not fully representative of conditions found in the eastern Middle Atlas and in the High Atlas, where unfavorable climatic conditions apparently threaten cedar populations (M'Hirit 2006).
Our research nevertheless gives insight into the ways combined pastoral and woodcutting activities affect Mediterranean forest structure in interaction with environmental parameters and natural ecosystem dynamics. However, more research is still needed for improved management. If C. atlantica forests are not likely to disappear in the next decades, research should aim at understanding and quantifying past changes and inferring their future. The understanding of recent changes could be completed, relying on a GIS diachronic analysis of old photographs, surveys with local actors, or new dendroecological or palynological records. Research could also contribute to understand how environmental factors interact with human activities to pinpoint the mechanisms responsible for unexplained differences and to assess perturbation regimes best for forest dynamics and biodiversity (Plieninger et al. 2003). More work is also needed in other C. atlantica stands and other Mediterranean forests. In particular, research could look at the water potential's role on height limitation during droughts (Koch et al. 2004) and the role of competition. In the Mediterranean region, ongoing research concerning the ecological importance and impacts of human activities will provide the technical understanding of sustainable ecosystem management. But environmental problems are societal problems, and management will never be effective if such research does not team up with socio-geographical studies to help understand why in ecosystems such as cedar forests conservation policies have so far mostly failed.