Understanding patterns and scales driving stand structure and diversity is a main topic of forest ecology (Legendre & Legendre, 1998). Forest dynamics drives the trees´ spatial patterns, while this spatial structure also determines to a large degree the properties of the whole forest ecosystem (Gadow et al., 2012). Specific processes generate particular patterns of growth and regenerationleading to upcontrasting structural arrangements (Gadow et al., 2012). Spatial pattern influences key aspects of ecosystem functions, including disturbance resistance and recovery dynamics, regeneration, snow retention, and habitat quality, as has been studied in fire-prone pine and mixed-conifer forests (Churchill et al., 2013). Hence, spatially explicit approaches are widely recommended for ecosystem management (Gadow et al., 2012). In the case of endangered tree species, spatial pattern analyses mightbe useful to guide conservation purposes (Raventós et al., 2011).

Consequently, stand structure and diversity has received increasing interest from forest managers in the recent decade (Schall et al., 2018). Multiple stand structural attributes characterize stand types, affecting management and driving the stand productivity (Schall et al., 2018). Fine-scale spatial structure is generally assessed by point pattern analyses based on second-order statistics of the overall point-to-point distances in a mapped region (Legendre & Legendre, 1998). These statistics provide a key tool for identifying spatial pattern types and scales and allow concluding if a pattern is clumped, regular or random (Legendre & Legendre, 1998; Wiegand & Moloney, 2004). Processes controlling spatial patterns of forest species has been investigated by point pattern analyses in temperate (Getzin et al., 2006, Carrer et al., 2018), tropical (Wiegand et al., 2007) and Mediterranean regions (Raventós et al., 2010, García-Cervigón et al., 2017). Here we focused on the Moroccan fir, Abies marocana Trab., a North African Tertiary relict species of the group of circum-mediterranean firs (Linares, 2011). This species has sometimes been considered as a subspecies of the Spanish fir, Abies pinsapo Boiss. Other authors suggested to split it into two species, Abies marocana and A. tazaotana Côzar ex Hug. Villar (Arista & Talavera, 1994; Farjon, 2010). More recent molecular and biometric studies have converged towards the separation of the Moroccan and Spanish firs (Sękiewicz et al., 2013, Dering et al., 2014). In the present study we adopted the classification of Fennane et al., (1999) considering Moroccan fir within a single endemic species.

Despite MF is included in a protected area, increasing human pressure, particularly by forest clearance for the expansion of Cannabis cultivation, logging, pollarding and grazing is still present. Furthermore, it is a drought-sensitive species, considered highly vulnerable to expected global warming during the second half of the 21st century (Sánchez-Salguero et al., 2017). Three decades ago, Melhaoui (1990) reported a lack of natural regeneration of MF and suggested two main causes: summer water stress and overgrazing by goats. Recently, Linares et al., (2011) reported a high regeneration dynamic of this species. Castello et al., (2016) found that MF population was structurally sustainable, whereas the combined fir populations in Spain and Morocco were structurally unsustainable due to excessive mortality. Actually, the species dynamic is not well understood while there is an urgent need to develop studies that inform conservation management.

The aim of this study was to assess spatial patterns of the Moroccan fir in forest stands in order to explore underlying processes. Given limited seed dispersal by wind in this fir and its tolerance to shade, we hypothesized intraspecific spatial aggregation patterns. These patterns of aggregation were tested by applying a univariate analysisand considering the location of stems (analysis 1). Spatial distribution of juveniles, relative to adult trees, was explored by a bivariate spatial analysis (analysis 2). Continuous size data (dbh) was included in marked point pattern analysis (analysis 3) to test the effect of proximity assuming that the dbh of two neighboring points tend to be smaller if they exhibit a competitive interaction (negative correlation), while they tend to have larger dbh than expected when a facilitative interaction is operating (positive correlation) (Getzin et al., 2006, 2008; Illian et al., 2008; Wiegand & Moloney, 2014). In the case of denser stands intraspecific competition may be strong and result in tree suppression or even mortality. We evaluated density-dependent suppression by random labeling analysis (analysis 4) of the relative distribution of dominant vs suppressed trees (Gavrikov & Stoyan, 1995). If suppressed trees are close to dominant ones we expect a one-sided intraspecific competition, whereas a symmetric competition would be inferred in the case of aggregation of suppressed individuals and their segregation from dominant trees (Raventós et al., 2010; Wiegand & Moloney, 2014).

Species and study site

Moroccan fir (MF) is located in a quaternary refuge with exceptional high moisture conditions, in the core area of Talassemtane National Park (TNP) in Northern Morocco. The area is included, since 2006, in the Mediterranean Intercontinental Biosphere Reserve of Spain and Morocco, established by the UNESCO’s Man and the Biosphere Program. Its forests have been recognized as a biodiversity hotspots and a centre of endemicity (Médail & Diadema, 2009). Moroccan fir is on the IUCN Red List as a globally threatened species in the category "Endangered" (Alaoui et al., 2011).

The study was carried out in the main distribution area of A. marocana at the TNP located in the western Rif Mountains, Chefchaoun, Morocco (Fig. 1). MF is distributed within two main forests: the Jbel Tazaout, accounting for about 300 ha and the Talassemtane range, accounting for about 3760 ha. MF is restricted to limestone and dolomitic substrates, at elevations between 1500 and 2000 m, on north, east and west slopes, in Mediterranean humid and perhumid bioclimates. In the upper limit of its area, MF is progressively replaced by the Atlas cedar Cedrus atlantica (Endl.) Carrière, while the Maghreb maritime pine Pinus pinaster Aiton subsp. hamiltonii var. maghrebiana Huguet del Villar dominates the lower part with the holm oak Quercus ilex L. and black pine Pinus nigra Arnold subsp. mauretanica (Maire and Peyer.) Heywood (Aafi, 2000). Soils are predominantly Alfisols and Mollisols developed on calcareous and dolomiticparent material (Benjelloun, 1993). Annual rainfall varies between 500 mm in the eastern valleys and more than 1500 mm on the high mountain peaks (Jbel Lakrâa 2156 m). The snow is present but occasional (Benabid, 2000). Mean temperature is around of 12-14 °C and extreme values for the mean of absolute maxima and absolute minima are 33°C and -3°C, respectively (Ghallab & Taïqui, 2015).

Field sampling

We sampled eight plots located in MF stands at the TNP with radius of 15 to 20 m depending on topography, for instance rock outcrops. In each plot, we recorded the distance of all trees to the center of plots and their azimuth, in order to determine their coordinates (x, y). A total of 908 trees (including 670 MF) were registered and measured, with a mean of 113 trees per plot. For each tree, we measured the diameter at breast height (dbh) of all individuals with dbh ≥ 2 cm. The status of trees was labeled as dominant individuals, suppressed individuals (trees green canopy only in the upper third), snags (standing dead trees), logs (downed dead wood) and stumps (cut-surface). We recorded the coordinates and topographic characteristics of each plot (elevation, aspect, and slope) in the field by means of a global positioning system (GPS). Observed signs of disturbance or human impacts were registered (fire, cutting, grazing, trampling). Thef ield work was carried out during the spring of 2010 and the summer of 2017.

Stand characterization

Density of stems (trees ha-1) with dbh ≥ 2 cm was calculated for MF in each of the eight plots as well as its basal area (m2 ha−1), mean dbh (cm) and the percentage of adults (dbh > 10 cm) and juveniles (2 ≤ dbh ≤ 10 cm). For other conifers and broad-leaved species, only stem densities are reported in this article. The analysis of these structural elements allows making an overall description of sampled plots. For each plot, the distribution of dbh was assembled on 5-cm classes. , was calculated both for woody species (pi corresponding to the proportion of trees being to the i-th species and S the number of woody species), and for dbh (pi representing the fraction of individuals belonging to the i-th dbh class and S the number of size classes) (Buongiorno et al., 1994).

Spatial pattern analyses

Analysis 1: We used the univariate pair-correlation function g(r) to measure both the intensity and scale thetree spatial pattern (Ripley, 1977; Stoyan & Stoyan, 1994; Wiegand & Moloney, 2004; Law et al., 2009). The g(r) is the expected number of points of the pattern at distance r from an arbitrary point of the pattern, divided by the pattern intensity λ in the plot (Wiegand & Moloney, 2004). It is calculated in the same way as Ripley’s K-function but with replacement of circles by rings, which gives to g(r) a non-cumulative character (Stoyan & Stoyan, 1994). In univariate analyses the Complete Spatial Randomness (CSR) is by far the most used null model (Velázquez et al., 2016) under which: g(r) = 1, g(r) > 1 for aggregation while g(r) < 1 for segregation. Because of the heterogeneity present in the plots, a Heterogeneous Poisson (HP) null model was used which assumes that the pattern intensity λ varies slightly with location (x, y) (Stoyan & Stoyan, 1994; Wiegand & Moloney, 2014). HP takes into account forest heterogeneity and thus reflects realistically tree-to-tree interactions (Carrer et al., 2018).

Analysis 2: We used a bivariate version of the pair-correlation function g12(r) to assess the relationship between adults (dbh > 10 cm) and MF juveniles (2 ≤ dbh ≤ 10 cm). It analyzes the association between all pairs of points of two patterns. The g12(r) is the density of pattern 2 points (juveniles in our case) at distance r from a typical point of pattern 1 (adults) divided by the density λ2 of pattern 2 points (Raventós et al., 2011). For this analysis, we tested the antecedent condition hypothesis (Wiegand & Moloney, 2004) by randomizing the position of the juveniles (pattern 2) and keeping the position of the adults (pattern 1) fixed. This null model is largely used in young-adult associations (Carrer et al. 2013, Velázquez et al. 2016) since it is applied for situations when two types of points are not created at the same time, but in sequence (Szmyt, 2014), which is the case of juvenile-adult relationships. Under antecedent condition, there is independence, attraction or repulsion for g12(r) = 1, g12(r) > 1 or g12(r) <1, respectively.

Analysis 3: In addition to point pattern locations, other information can be integrated in the analysis, usually called as “marks” such as tree size (dbh in our case) (Wiegand & Moloney, 2004). Effect of intraspecif ic interaction on individual growth was analyzed using univariate mark correlation function kmm(r) which quantifies the similitude or dissimilitude between the marks mi and mj of two neighboring fir trees i and j, separated by distance r (Illian et al., 2008; Wiegand & Moloney, 2014). It corresponds to the test function of size product f(mi, mj) = mi*mj normalized by division by the squared mean mark for all marks in the plot (Getzin et al., 2008). If kmm(r) > 1, the two trees tend to have larger dbh than average marks when they are nearby (e.g. positive correlation), while a value of kmm (r) < 1 indicates individuals tend to have smaller dbh (e.g. negative correlation), and if kmm(r) = 1 the two marks are independently assignedto individuals (Illian et al., 2008; Wiegand & Moloney, 2014). We used the independent marking null model that randomly shuffles the mark among all MF individuals (Wiegand & Moloney, 2014). It assumes that dbh values are independent of individual location (Stoyan & Stoyan, 1994).

Analysis 4: We used random labeling null-model to test partial-pair correlation function and mark connection function of two categories of MF trees: suppressed trees vs trees in good condition (hereafter referred to as dominant trees). As suggested by Goreaud & Pélissier (2003), random labeling null model is appropriate when such categories are considered to result mainly from different growing conditions affecting a posteriori a single population. The bivariate pair correlation has been defined in analysis 2 and it is applied here to test attraction vs repulsion among these two patterns (suppressed and dominant trees) under random labeling null model. The mark connection function p12(r) describes the probability of suppression of MF trees in dependence of the distance r to conspecific trees (Gavrikov & Stoyan, 1995; Illian et al., 2008). Under random labeling, p12(r) = p1 p2, where p2 is the proportion of suppressed trees among all MF individual. For p12(r) < p1 p2, suppressed trees and trees dominant trees are segregated at distance r within all individuals while larger values of p12 indicate a positive association (Gavrikov & Stoyan, 1995; Raventós et al., 2010).

For each type of analysis, data from the eight sampled plots were combined inside one average function (Illian et al., 2008; Wiegand & Moloney, 2014) to serve as statistical replicates. The eight sampled plots have been selected according to conditions of replicate analysis (same study area, within a limited range of environmental variation) from a total of 19 original plots. The common method is to estimate the statistical functions separately for each plot and then to aggregate the estimates in one mean overall test statistic (Illian et al. 2008). Practically, we performed data analysis for each plot separately, and then we combined the results into a single summary statistic and produce the associated simulation envelopes. Programita ® (htpps://programita.org), a software for spatial point pattern analysis in Ecology, uses an aggregation formula in which the estimator of the O-ring statistic is combined with that of the intensity (Equations 3.107 and 3.108 in Wiegand & Moloney, 2014, p. 249-250) to obtain the final aggregated estimator of each test function.

For all analyses and for each null model we performed 199 Monte Carlo simulations of the envelope limits between fifth lowest and highest values (Wiegand & Moloney, 2014). If the curve of observed function fall outside the upper simulation envelope there is a positive association (aggregation or attraction), while if the empirical curve is below the lower simulation envelope a negative association is expected (segregation or repulsion). Finally,the patterns are independent if the observed function is fitting inside simulation envelopes. To assess the ability of the null models to reflect the data we used a goodness-of-fit test GoF (Loosmore & Ford, 2006).All spatial pattern analyses were done with Programi-ta software (Wiegand & Moloney, 2014).

Stand composition and characteristics

Moroccan fir density, basal area and tree size and other structural characteristics across the eight studied plots aregiven in Table 1. MF Maximum density and highest basal area were recorded in plot 3 with 1386 trees ha-1 and 68.9 m2 ha-1, respectively; the minimum ones in plot 7 and plot 6 with 392.6 trees ha-1 and 19.1 m2 ha-1, respectively. The percentage of juveniles was higher than that of adults in plots situated at higher altitudes. Mean dbh of MF was the highest in the low elevation plot 2 (20.6 cm) and smaller in plots 3 to 6 (similar mean values around 13 cm) (Table 1).

In terms of MF status, the densest plot (plot 3) contained the highest density of suppressed and dead trees (580.0 and 155.6 trees ha-1, respectively, Table 1). Conversely, less dense plots generally showed lower density values of dead wood which did not exceed 8.0 trees ha-1in plot 1 and was absent in plot 6. Plot 8 contained the highest density of stumps (95.5 stem ha-1) followed by plots 7 and 5 (67.2 and 63.7 trees ha-1, respectively, Table 1).

There were other woody species that were classified into two groups (Table 1): conifers (mainly Cedrus atlantica, Juniperus oxycedrus and Pinus. mauretanica) and broad-leaved species (mainly Acer granatense and Crataegus laciniata). Cedrus atlantica was the second dominant species (after MF). Its density reached 396.1 trees ha-1 in plot 3. Despite its low density, Acer granatense was found in all except one plot. Other important species occurred to a lesser extent, such as Pinus nigra, Juniperus oxycedrus and Crataegus laciniata. The plot 6 comprised almost all these species and thus presented the highest Shannon’s species diversity index (H' species, Table 1), while plots 7 and 8 had the lower values of H' species (with 0.25 and 0.5, respectively). The highest value of H’ (dbh) was found in plot 2 (H'dbh=2.19), whereas the smallest values were recorded in plots 8 and 6 (1.68 and 1.69, respectively). Diameter distribution of MF in all sampled plots showed reverse-J shape (Fig. 2) with a high number of trees in the smaller size classes but wider size classes were relatively rare. All plots showed diameter distribution peak in the smallest dbh class (2–5 cm), and all had stems until the diameter range of 30–35 cm. Plots 1, 2 and 3 presented the larger diameter distribution, while plot 4 did not contain individuals larger than 55 cm of dbh.

Spatial pattern of Moroccan fir

The number of Moroccan fir individuals included in spatial pattern analyses is reported between brackets for each plot in Table 1 and the locations of all measured trees by plot are reported in Figure S1 [suppl.]. Result of univariate pair correlation function using all MF trees is shown in Figure 3a. The g(r) of the replicated plots indicated asignificant aggregation (p = 0.005) up to 3 m. Individuals of MF have more conspecifics in their close neighborhood than expected by the null model. Bivariate analysis of the relationship between MF juveniles (2 ≤ dbh ≤ 10 cm) and adults (dbh > 10 cm) showed a strong significant attraction between the two size classes in the combined plots up to 2 m (p = 0.005, Fig. 3b). Juveniles were clearly more aggregated around adult trees than expected by the null model. The mark-correlation functions kmm(r) showed a non-significant negative correlation (p = 0.015) up to a distance of 2 m and at 9-11 m. This indicated that sizes of MF individuals were independently distributed from their location (Fig. 4a). Finally, for the random labeling analysis, the pair correlation function g12(r) showed a weak but significant departure from the null model at 2.5-3.5 m (p = 0.005, Fig. 4b). This indicates that there is repulsion between suppressed and dominant individuals. The mark connection function p12(r) revealed similar results; a negative association between dominant and suppressed trees is shown at a small-scale of 2.5–3.5 m (p = 0.005, Fig. 4c). Consequently, when randomly selecting two points with this scale r there was a high probability that both were dominant or that both were suppressed.

The density of young trees was higher than that of adults in all except one plot and reached 961.9 trees ha-1 in plot 3 (Table 1). The frequency distribution of tree dbh classes showed an overall J-reverse shape reflecting a good regeneration capacity with a prompt decline in the abundance of larger-size trees (Fig. 2). Similar results were found by Linares et al., (2011) and recently by Navarro-Cerrillo et al. (2020) in TNP. Tree density and basal area were higher at high elevation and were associated with increasing rates of mortality. Besides, the percentageof young trees increased as elevation did, while mean tree diameter was relatively higher at lower elevation (Table 1). A similar decrease of tree height as elevation rises has been reported by Benabid (1982). The resulting diversity of tree sizes is lower at low altitude. Regarding species richness, overall diversity of woody species was relatively higher at middle and lower elevations (plots 4 and 6).

MF attained relatively higher mean dbh values in less dense stands (plots 2 and 7) which have lower values of Shannon’s dbh index (H’ dbh, Table 1). Similar relationships between size diversity, wood species diversity and stand density were found by Linares et al., (2011). Compared to Spanish fir forests, MF forests are characterized by higher species diversity (Aafi, 2000; Benabid, 2000; Linares et al., 2011) due to the presence of several tree species mainly C. atlantica, A. granatense, P. nigra and J. oxycedrus. The presence of these species with their structural features descripted in the present study was similar to the results found by Navarro-Cerrillo et al. (2020).

Nevertheless, our results also support other environmental factors affecting MF forest structure and diversity.Hence, stumps frequency, especially of small sized trees indicating human perturbation, mainly at low elevation (plots 8, 7, 6, 5 and 2). In the Rif, forests have deeply suffered from an abusive exploitation since the protectorate period (Benabid, 2002) which was reflected by the absence of very large size classes (Fig. 2). However, before this period, the exploitation by local populations was relatively less important although it has considerably increased over last decades (Benabid, 2002). In a comparison between Moroccan and Spanish fir stands, Linares et al., (2011) reported a widespread logging activity by local inhabitants in Talassemtane forest. Furthermore, anthropogenic pressure, especially the extension of Cannabis crops, remains a major threat to the maintenance of MF forests (Benabid, 2000; Taïqui, 2005). Previous studies reported the predominance of these crops at low altitudes; however, at present, Cannabis crops are increasingly becoming the main driver of deforestation of MF and other mountainous forests in Northern Morocco.

We found a significant spatial aggregation of Moroccan fir trees, at a fine scale of up to 3 m, as showed by theunivariate pair-correlation function (Fig. 3a). This result is refined by the bivariate analysis of two size classes; MF young individuals are mainly concentrated around adult trees at small-scale (up to 2 m; Fig. 3b). Additionally, mar-ked point pattern analyses showed that individual tree dbh is independently distributed from location, which indicate that aggregation does not affect tree size (Fig. 4a). As suggested for other Mediterranean ecosystems (De Luis et al., 2008), the observed MF pattern may be produced by a combination of various mechanisms including its limited seed dispersal, shade tolerance and the ecological effects of spatial heterogeneity of resources and conditions. The general context of climatic stress and human pressure should be also considered.

Spatial patterns of tree species usually show a scale-dependent aggregation that is generated by seed dispersal limitation (De Luis et al., 2008; Zhang et al., 2013). Wind dispersal can cause short distance dissemination away from source leading to seed shadow around adult trees (Wiegand et al., 2007). Seeds that did not dispersed far from adult trees are accumulated under the crowns so that few juveniles are found far from the adult stems (LeMay et al., 2009; Raventós et al., 2011; García-Cervigón et al., 2017). Like other Mediterranean firs and pines (García-Cervigón et al., 2017; Abellanas & Pérez Moreno, 2018), MF seeds are wind-dispersed. Melhaoui (1990) estimated that the dispersal distance of MF is almost always less than 5 m. Our finding supports this estimation since the aggregation of juveniles around MF adults occurred up to 2 m (Fig. 3b).

However, seedling establishment at short distance from adults depends on species preferences and it is usually thought to be affected by spatial heterogeneity and its effects on habitat quality (Zhang et al., 2013). For example, shade tolerance or intolerance abilities can play an important role in species spatial pattern (Abellanas & Pérez-Moreno, 2018). But unlike wet sites where seedlings development is strongly influenced by light (Chave, 2000), the limiting factor at dry sites seems to be moisture where large trees offer moister microsites to small trees (LeMay et al., 2009). In our case, like the majority of fir species, MF is a shade-tolerant species that regenerate under canopy layer (Baumer, 1977). Previous studies on MF also found that the moisture factor plays a crucial role in natural regeneration (Melhaoui, 1990). Clumping of MF juveniles under adult trees is due to the favorable moisture conditions below canopy and shade-tolerance ability of juveniles. Taking separately, almost all plots showed aggregated pattern of MF trees. The exceptions were plots 6 and 8 for the univariate analysis and plots 5 and 8 for the bivariate one, which showed random and independent patterns, respectively (Fig. S2 and S3 [suppl.]). These departures from the general findings might be explained by a higher anthropogenic pressure reflected by the high amount of stumps in these plots. Indeed, many studies reported that anthropogenic activities such as logging can modify the spatial structure of the stand (Motta & Edouard, 2005).

In addition to seed dispersal limitation and environmental heterogeneity, biotic interactions play an important role in the dynamics of an ecosystem and can induce specific structures (Goreaud & Pélissier, 2003). Previous studies have used second order statistics to study the relationships between spatial patterns and biotic processes, and have reported positive, independent and negative correlations between tree positions and sizes. In a study of Pistacia atlantica in a semi-arid region of Iran (Erfanifard et al., 2018), bivariate analysis revealed a spatial aggregation of saplings around females while mark analysis using kmm(r) showed a negative correlation between dbh and position. In this case, aggregation negatively affects size. In a large study area covering 30 ha in a temperate old-growth forest in China, Zhang et al., (2013) identif ied mostly independence of tree dbh at greater distances using kmm(r), but highly significant negative correlation at short distances for the majority of species and some cases of significant short-distance attraction of dbh associated with the light-demanding shrub species. Jácome-Flores et al., (2016) found at Martinazo site in Spain a clustered distribution of dwarf palms Chamaerops humilis individuals while their size did not show a spatial correlation. They suggested that this result can be seen when using small sample sizes, which produce noticeably wider simulation envelopes. Getzin et al., (2008) showed no correlation of the studied marks (dbh, crown area and upper crown area) using kmm(r) at all scales in a low-density deciduous stand in Germany which imply low overall competition effect. Similar results were found for Pinus syivestris in Poland where no correlation was found between trees either using diameter or height of pine individuals (Erfanifard & Stereńczak, 2017). Law et al., (2009) reported that the absence of correlation may suggest that the size is much more equally distributed over space than the individuals are themselves. In our study, mark correlation analysis of Abies marocana trees in Talassemtane forest showed a non-significant negative correlation of tree sizes up to 2 m, as well as at 9-11 m where a negligible negative correlation appeared (Fig. 4a). This non-significant negative correlation may be attributed to the occurrence of a weak competitive effect at small distance. Negative correlation indicates aggregation of young individuals and may be related to clumped dispersal and establishment patterns, which are absent at later developmental stages due to competition and mortality.

At the plot level, we found no correlation of tree sizes except for plot 2 which showed a very weak and non-significant negative correlation up to a very short distance (< 0.5 m) (Fig. S4 [suppl.]). The significant pattern detected by the average function at small scale based on replicated pattern could be reinforced by the additional effect of weak negative correlation tendency observed in plots 5, 7 and 8 (where the kmm(r) expected value was closed to the lower simulation envelope). However, the facilitative effect recorded among young and aged trees may hide the competitive effect occurring within groups of regenerating individuals.

Furthermore, low density of suppressed trees as well as their spatial pattern confirms weak competitive interaction and self-thinning process. In fact, the result of random labeling either by g12(r) or p12(r) showed a weak but significant repulsion between suppressed and dominant individuals at 2.5-3.5 m (p = 0.005, Fig. 4b and 4c). The highest density of suppressed trees in association with a higherrate of dead trees in the densest plot 3 (Table 1) may indicates a density-dependent suppression process but which could be significantly detected by test functions only at a limited scale. This negative association appears to be a small-scale pattern only emerging within a 2.5–3.5 m range which reflects a fine patchy pattern where small and almost distinct groups of dominant and suppressed trees occur separately within the stands (Carrer et al., 2013). As most of the suppressed trees were small trees, and MF regeneration was abundant in most plots, the spatial negative correlation in tree size could contribute sporadically to suppression of tree growing close to each other in denser stands. Accordingly, MF juveniles clumped in the shelter of the adults and beneficiate from a positive recruitment interaction which seems to be more important than an intra-specific competitive effect.

Here, the random labeling considered for each plot showed the distinct behavior of the four last plots (Fig. S5 [suppl.]), and thus support the previous observations in the three other analysis types (Fig. S2, S3 and S4 [suppl.]). Indeed, a negative association between dominant and suppressed trees is found in plots 5, 8 and particularly 7 which extends to 4.5 m, in addition to a significant positive association at distance of > 14 m shown in this last plot. At the small scale, suppressed trees appear to be separated from dominant ones while they are more grouped at larger scale.

Under current climate conditions and human pressure in TNP (Taïqui, 2005; Linares et al., 2011; Ghallab & Taïqui, 2015), seedlings and saplings of Moroccan fir regenerate and maintain themselves in favorable moisture conditions under the canopy of adult trees thanks to their shade tolerance and the absence of significant intraspecific competitive effect. Spatial intraspecific negative interaction may sometimes be undetectable because MF can withstand competition for light since it is a shade tolerant species (Baumer, 1977; Abellanas & Pérez-Moreno, 2018). This suggests that intraspecific competition did not appear to be a major process that negatively affects MF growth in the studied stands, which is just the opposite of what was found for the Spanish fir counterpart (Linares et al., 2010). In contrast, the facilitative interaction where parents provide shelter from abiotic stress and protection against biotic consumption was more relevant. Although our study has provided insights into the spatial patterns of Moroccan fir, there are some limitations regarding the relative number of points used within each plot, particularly in the case of analyses involving two pattern types. Indeed, Wiegand & Moloney (2014) recommended the use of 50-70 individuals as a minimum in spatial point analyses. However, Velázquez et al. (2016) reported that many studies used relatively few individuals in point pattern analysis (< 100). For instance, Camarero et al. (2000) used at least 15 individuals, Zenner & Peck (2009) used 20 individuals. More recently, Janík et al. (2016) and Cordero et al. (2016) used much less trees (10 and 4, respectively). In our case, we had a sufficient number of trees in univariate analyses, and despite that the available number of trees was split into two classes in the bivariate and random labeling analyses, this limitation was minimized by using replicate analysis with data from the eight sampled plots that were combined into one average function (Illian et al., 2008; Wiegand & Moloney, 2014). This method is a straight-forward extension of the typical analysis of separate plots (De Luis et al., 2008) and allows increasing the sample size and the statistical power (Illian et al., 2008). Although being considered as an extension of the common analysis of isolated plots (De Luis et al., 2008), replicate analysis was rarely used by ecologists in spatial point patterns (but see Raventós et al., 2011).

The results of our explanatory study allowed us to describe the overall spatial pattern of the Moroccan endemic f ir. Analysis of the main MF structural characteristics in the sampled plots indicated a remarkable regeneration dynamic. However, the anthropic pressure remains a serious problem which constitutes a constraint to MF community preservation. Moroccan fir showed strong positive intraspecific association and juveniles were significantly attracted to adult trees. We highlight the importance of shade tolerance, seed dispersal and microsite heterogeneity in the observed facilitative intraspecific effect. However, more investigation is needed to test subsequent hypothesis, in particular by integrating other species (e.g. Cedrus atlantica). Understanding spatial structure and dynamics of trees, especially endemic species, is critical particularly in the context of global change.

Aafi A, 2000. Floristic diversity of Morocco's fir ecosystem (Abies maroccana Trab.) (Talassemtane National Park). Nature et Faune 18 (1): 15-19.Abellanas B, Pérez-Moreno PJ, 2018. Assessing spatial dynamics of a Pinus nigra subsp. salzmannii natural stand combining point and polygon patterns analysis. For Ecol Manage 424: 136-153.https://doi.org/10.1016/j.foreco.2018.04.050Alaoui ML, Knees S, Gardner M, 2011. Abies pinsapo var. marocana, Moroccan Fir. The IUCN Red List of Threatened Species 2011: E.T34126A9841418. URL:http://www.iucnredlist.org/details/34126/0Arista M, Talavera S, 1994. Phenology and anatomy of the reproductive phase of Abies pinsapo Boiss. (Pinaceae). B J Linn Soc 116: 223-234.https://doi.org/10.1111/j.1095-8339.1994.tb00431.xBaumer M, 1977. Le sapin du Maroc. R.F .F ., XXIX: 343-354.https://doi.org/10.4267/2042/21153Benabid A, 1982. Etude phytoécologique, biogéographique et dynamique des associations sylvatiques du Rif occidental. Problémes posés par la réforestation et l'aménagement des peuplements forestiers actuels. Doctoral thesis, Université St Jérome-Marseille.Benabid A, 2000. Flore et écosystèmes du Maroc: Évaluation et préservation. Ibis Press, Paris, France. 359 pp.Benabid A, 2002. Le Rif et le Moyen-Atlas (Maroc): Biodiversité, Menaces, Préservation. Mountains High Summit Conference for Africa. UNEP, Kenya (Nairobi) May 6-10.Benjelloun H, 1993. Soil genesis, classification, and nitrogen cycling in forest ecosystems of the Northwestern Rif region of Morocco. Doctoral thesis, Oregon State University, USA.Buongiorno J., Dahir S., Lu H.C., and Lin C.R., 1994. Tree size diversity and economic returns in uneven-aged forest stands. For Sci 40: 83-103.Camarero JJ, Gutiérrez E, Fortin MJ, 2000. Spatial pattern of subalpine forest-alpine grassland ecotones in the Spanish Central Pyrenees. For Ecol Manage 134: 1-16.https://doi.org/10.1016/S0378-1127(99)00241-8Carrer M, Soraruf L, Lingua E, 2013. Convergent space - time tree regeneration patterns along an elevation gradient at high altitude in the Alps. For Ecol Manage 304: 1-9.https://doi.org/10.1016/j.foreco.2013.04.025Carrer M, Castagneri D, Popa I, Pividori M, Lingua E, 2018. Tree spatial patterns and stand attributes in temperate forests: The importance of plot size, sampling design, and null model. For Ecol Manage 407: 125-134.https://doi.org/10.1016/j.foreco.2017.10.041Castello J D, Cale JA, Angelo CMD, Linares JC, 2016. Baseline Mortality Analysis Reveals Legacy of Contrasting Land Use Practices on the Structural Sustainability of Endangered Moroccan and Spanish Mountain Forests. Forests 7: 1-12.https://doi.org/10.3390/f7080172Chave J, 2000. Dynamique spatio-temporelle de la forêt tropicale. Ann Phys Fr 25: 1-184.https://doi.org/10.1051/anphys:200006001Churchill DJ, Larson AJ, Dahlgreen MC, Franklin JF, Hessburg PF, Lutz JA, 2013. Restoring forest resilience: From reference spatial patterns to silvicultural prescriptions and monitoring. For Ecol Manage 291: 442-457.https://doi.org/10.1016/j.foreco.2012.11.007Cordero I, Jiménez MD, Delgado JA, Villegas L, Balaguer L, 2016. Spatial and demographic structure of tara stands (Caesalpinia spinosa) in Peru: Influence of present and past forest management. For Ecol Manage 377: 71-82.https://doi.org/10.1016/j.foreco.2016.06.034Dering M, Sękiewicz K, Boratyńska K, Litkowiec M, Iszkuło G, Romo A, Boratyńska A, 2014. Genetic diversity and inter-specific relations of western Mediterranean relic Abies taxa as compared to the Iberian A. alba. Flora 209: 367-374.https://doi.org/10.1016/j.flora.2014.03.011De Luis M, Raventós J, Wiegand T, Gonzailez-hidalgo JC, 2008. Temporal and Spatial Differentiation in Seedling Emergence May Promote Species Coexistence in Mediterranean Fire-Prone Ecosystems. Ecography 31: 620-629.https://doi.org/10.1111/j.0906-7590.2008.05433.xErfanifard Y, Stereńczak K, 2017. Intra- and interspecific interactions of Scots pine and European beech in mixed secondary forests. Acta Oecol 78: 15-25https://doi.org/10.1016/j.actao.2016.12.002Erfanifard Y, Nguyen HH, Schmidt JP, Rayburn A, 2018. Fine-scale intraspecific interactions and environmental heterogeneity drive the spatial structure in old-growth stands of a dioecious tree. For Ecol Manage 425: 92-99.https://doi.org/10.1016/j.foreco.2018.05.041Farjon A, 2010. A Handbook of the World's Conifers. Koninklijkr Brill NV, Leiden.https://doi.org/10.1163/9789047430629Fennane M, Ibn Tattou M, Mathez J, Ouyahya A, El Oualidi J, 1999. Flore pratique du Maroc. Pteridophyta, Gymnospermae, Angiospermae (Lauraceae - Neuradaceae). Travaux de l'Institut Scientifique série Botanique, Agdal, Rabat. 559 pp.Gadow K, Zhang CY, Wehenkel C, Pommerening A, Corral-Rivas J, Korol M, Myklush S, Hui GY, Kiviste A, Zhao XH, 2012. Forest Structure and Diversity. In: Continuous Cover Forestry; Pukkala T, Gadow VK, (Eds). Springer Netherlands: Dordrecht, The Netherlands pp 29-83.https://doi.org/10.1007/978-94-007-2202-6_2García-Cervigón AI, Velázquez E, Wiegand T, Escudero A, Olano JM, 2017. Colonization in Mediterranean old-fields: the role of dispersal and tree - tree interactions. J Veget Sci, 28, 627-638.https://doi.org/10.1111/jvs.12500Gavrikov V, Stoyan D, 1995. The use of marked point processes in ecological and environmental forest studies. Environ Ecol Stat 2: 331-344.https://doi.org/10.1007/BF00569362Getzin S, Dean C, He F, Trofymow JA, Wiegand K, Wiegand T, 2006. Spatial patterns and competition of tree species in a Douglas-fir chronosequence on Vancouver Island. Ecography 29: 671-682.https://doi.org/10.1111/j.2006.0906-7590.04675.xGetzin S, Wiegand K, Schumacher J, Gougeon A, 2008. Scale-dependent competition at the stand level assessed from crown areas. For Ecol Manage 255: 2478-2485.https://doi.org/10.1016/j.foreco.2008.01.007Ghallab A, Taïqui L, 2015. Modélisation de la distribution spatiale des paramètres bioclimatiques dans la région Tanger- Tétouan (Rif occidental). Eur Sci J 11: 1857-7881.Goreaud F, Pélissier R, 2003. Avoiding misinterpretation of biotic interactions with the intertype-function: population independence vs. random labeling hypothesis. J Veget Sci 14: 681-692.https://doi.org/10.1111/j.1654-1103.2003.tb02200.xIllian J, Penttinen A, Stoyan H, Stoyan D, 2008. Statistical Analysis and Modelling of Spatial Point Patterns. International Statistical Review. John Wiley, New York, USA. 534 pp.https://doi.org/10.1002/9780470725160Janík D, Král K, Adam D, Hort L, Samonil P, Unar P, Vrska T, Mcmahon S, 2016. Tree spatial patterns of Fagus sylvatica expansion over 37 years. For Ecol Manage 375: 134-145.https://doi.org/10.1016/j.foreco.2016.05.017Jácome-Flores ME, Delibes M, Wiegand T, Fedriani JM, 2016. Spatial patterns of an endemic Mediterranean palm recolonizing old fields. Ecol Evol 6: 8556-8568https://doi.org/10.1002/ece3.2504Law R, Illian J, Burslem DFRP, Gratzer G, Gunatilleke CVS, Gunatilleke IAUN, 2009. Ecoogical information from satial patterns of plants: Insights from point process theory. J Ecol 97: 616-628.https://doi.org/10.1111/j.1365-2745.2009.01510.xLegendre P, Legendre L, 1998. Numerical Ecology (2nd English edn). Amsterdam, The Netherlands: Elsevier. 853 pp.LeMay V, Pommerening A, Marshall P, 2009. Spatio-temporal structure of multi-storied, multi-aged interior Douglas fir (Pseudotsuga menziesii var. glauca) stands. J Ecol 97: 1062-1074.https://doi.org/10.1111/j.1365-2745.2009.01542.xLinares JC, 2011. Biogeography and evolution of Abies (Pinaceae) in the Mediterranean Basin: The roles of long-term climatic change and glacial refugia. J Biogeogr 38: 619-630.https://doi.org/10.1111/j.1365-2699.2010.02458.xLinares JC, Camarero JJ, Carreir JA, 2010. Competition modulates the adaptation capacity of forests to climatic stress: insights from recent growth decline and death in relict stands of the Mediterranean fir Abies pinsapo. J Ecol 98: 592-603.https://doi.org/10.1111/j.1365-2745.2010.01645.xLinares JC, Carreira JA, Ochoa V, 2011. Human impacts drive forest structure and diversity. Insights from Mediterranean mountain forest dominated by Abies pinsapo (Boiss.). European J For Res 130: 533-542.https://doi.org/10.1007/s10342-010-0441-9Loosmore NB, Ford ED, 2006. Statistical inference using the g or k point pattern spatial statistics. Ecology 87: 1925-1931.https://doi.org/10.1890/0012-9658(2006)87[1925:SIUTGO]2.0.CO;2Médail F, Diadema K, 2009. Glacial refugia influence plant diversity patterns in the Mediterranean Basin. J Biogeogr 36: 1333-1345.https://doi.org/10.1111/j.1365-2699.2008.02051.xMelhaoui Y, 1990. Etude phytoécologique, productivilé et classes de croissance du Sapin du Maroc (Abies maroccana Trab.). Problématique de la régénération naturelle des peuplements de la sapiniére marocaine. Doctoral thesis, Aix-Marseille III Univ, Aix-Marseille, France.Motta R, Edouard JL, 2005. Stand structure and dynamics in a mixed and multilayered forest in the Upper Susa Valley, Piedmont, Italy. Can J For Res 36: 21-36.https://doi.org/10.1139/x04-153Navarro-Cerrillo RM, Manzanedo RD, Rodriguez-Vallejo C, Gazol A, Palacios-Rodríguez G, Camarero JJ, 2020. Competition modulates the response of growth to climate in pure and mixed Abies pinsapo subsp. Maroccana forests in northern Morocco. For Ecol Manage 459: 117847.https://doi.org/10.1016/j.foreco.2019.117847Raventós J, Wiegand T, De Luis M, 2010. Evidence for the spatial segregation hypothesis : a test with nine-year survivorship data in a Mediterranean shrubland. Ecology: 91: 2110-2120.https://doi.org/10.1890/09-0385.1Raventós J, Mujica E, Wiegand T, Bonet A, 2011. Analyzing the Spatial Structure of Broughtonia cubensis (Orchidaceae) Populations in the Dry Forests of Guanahacabibes, Cuba. Biotropica 43: 173-182.https://doi.org/10.1111/j.1744-7429.2010.00684.xRipley BD, 1977. Modelling Spatial Patterns. Journal of the Royal Statistical Society 39: 172-212.https://doi.org/10.1111/j.2517-6161.1977.tb01615.xSánchez-Salguero R, Camarero JJ, Carrer M, Gutiérrez E, Alla AQ, Andreu-Hayles L, Hevia A, Koutavas A, Martínez-Sancho E, Nola P, et al., 2017. Climate extremes and predicted warming threaten Mediterranean Holocene firs forests refugia. Proceedings of the National Academy of Sciences of the United States of America 114: 10142-10150.https://doi.org/10.1073/pnas.1708109114Schall P, Schulze ED, Fischer M, Ayasse M, Ammer C, 2018. Relations between forest management, stand structure and productivity across different types of Central European forests. Basic App Ecol: 32: 39-52.https://doi.org/10.1016/j.baae.2018.02.007Sękiewicz K, Sękiewicz M, Jasińska AK, Boratyńska K, Iszkuło G, Romo A, Boratyński A, 2013. Morphological diversity and structure of West Mediterranean Abies species. Plant Biosyst 147: 125-134.https://doi.org/10.1080/11263504.2012.753130Stoyan D, Stoyan H, 1994. Fractals, random shapes and point fields. Methods of geometrical gtatistics. John Wiley & Sons, Chichester, UK. 394 pp.Szmyt J, 2014. Spatial statistics in ecological analysis: From indices to functions. Silva Fennica 48: 1-31.https://doi.org/10.14214/sf.1008Taïqui L, 2005. Evolution récente de la structure du paysage du Bassin de Chefchaouen (1958-1986). Mutations Des Milieux Ruraux Dans Les Montagnes Rifaines (Maroc), Série Etudes Spatiales 2: 1-19.Velázquez E, Martínez I, Getzin S, Moloney KA, Wiegand T, 2016. An evaluation of the state of spatial point pattern analysis in ecology. Ecography 39: 001-014.https://doi.org/10.1111/ecog.01579Wiegand T, Gunatilleke S, Gunatilleke N, 2007. Species Associations in a Heterogeneous Sri Lankan Dipterocarp Forest. Am Nat 170: E77-E95.https://doi.org/10.1086/521240Wiegand T, Moloney KA, 2014. Handbook of Spatial Point- Pattern Analysis in Ecology. Chapman and Hall/CRC Press, Boca Raton, FL. 212-234 pp.https://doi.org/10.1201/b16195Wiegand T, Moloney KA, 2004. Rings, circles, and null-models for point pattern analysis in ecology. Oikos 2: 209-229.https://doi.org/10.1111/j.0030-1299.2004.12497.xZenner EK, Peck JE, 2009. Characterizing structural conditions in mature managed red pine: Spatial dependency of metrics and adequacy of plot size. For Ecol Manage 257: 311-320.https://doi.org/10.1016/j.foreco.2008.09.006Zhang C, Wei Y, Zhao X, Gadow KV, 2013. Spatial characteristics of tree diameter distributions in a temperate old-growth forest. PLoS ONE 8: e58983.https://doi.org/10.1371/journal.pone.0058983