Characterizing spatial structure of Abies marocana forest through point pattern analysis

Aim of study: Understanding small-scale patterns caused by stochastic factors or community interactions driving forest structure and diversity of Moroccan fir Abies marocana Trab. Area of study: Talassemtane fir forest, Talassemtane National Park, Rif Mountains, northern Morocco. Material and methods: Eight plots representative of the structural variability of A. marocana forests were selected, and all tree individuals with diameter at breast height (dbh) ≥2 cm were mapped and measured. We performed four types of spatial point pattern analyses: (1) Univariate analyses to explore the overall trees spatial pattern, (2) bivariate analyses to assess the spatial relationship between juveniles and A. marocana adults, (3) correlation between tree sizes (dbh) and distance between points pairs using the univariate mark correlation function, and (4) random labeling analysis between dominant and suppressed Moroccan fir individuals to assess competition patterns. Main results: We found a strong spatial aggregation of fir individuals and a positive intraspecific association between juveniles and adult trees. However, there were weak but significant distance-dependent effect on tree size and density-dependent effect on suppression pattern. Research highlights: Shade-tolerance, seed dispersal and/or microsite heterogeneity might play important roles in the observed fir pat -terns. Our results provide a basic knowledge on within-stand Moroccan fir spatial distribution, with implications for adaptive management of these relic forests, and prompting to further research to test advanced hypotheses.


Introduction
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 regeneration leading 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 might be 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 analysis and 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(Getzin et al., , 2008Illian 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 dolomitic 3 Spatial structure of Abies marocana parent 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). The field 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 (m 2 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 the tree 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 g 12 (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 g 12 (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 g 12 (r) = 1, g 12 (r) > 1 or g 12 (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 intraspecific interaction on individual growth was analyzed using univariate mark correlation function k mm (r) which quantifies the similitude or dissimilitude between the marks m i and m j 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(m i , m j ) = m i *m j normalized by division by the squared mean mark for all marks in the plot (Getzin et al., 2008). If k mm (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 k mm (r) < 1 indicates individuals tend to have smaller dbh (e.g. negative correlation), and if k mm (r) = 1 the two marks are independently assigned to 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 p 12 (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, p 12 (r) = p 1 p 2 , where p 2 is the proportion of suppressed trees among all MF individual. For p 12 (r) < p 1 p 2 , suppressed trees and trees dominant trees are segregated at distance r within all individuals while larger values of p 12 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, Spatial structure of Abies marocana 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 goodnessof-fit test GoF (Loosmore & Ford, 2006).

Stand composition and characteristics
Moroccan fir density, basal area and tree size and other structural characteristics across the eight studied plots are given in Table 1. MF Maximum density and highest basal area were recorded in plot 3 with 1386 trees ha -1 and 68.9 m 2 ha -1 , respectively; the minimum ones in plot 7 and plot 6 with 392.6 trees ha -1 and 19.1 m 2 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 -1 In the present study, we considered as juveniles individuals with 2 ≤ dbh ≤ 10 cm and as adults those with dbh > 10 cm. For Moroccan fir, we estimated tree status as: dead individuals (snags and logs), suppressed and stumps. For the other species, only the total density is shown. Number of Moroccan fir individuals included in spatial analyses is shown between brackets for each plot in 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 a significant 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 k mm (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 g 12 (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 p 12 (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.

Discussion
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 percentage Spatial structure of Abies marocana of 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 the univariate 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, marked 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 su-ggested 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. Random labeling analysis of firs (living trees vs suppressed trees) using partial-pair correlation functions g 12 (r) and (c) mark connection function p 12 (r). Dotted line indicates the empirical curve, gray line indicates the expected value under the null model (independent marking for k mm (r) and random labeling for g 12 (r) and p 12 (r)) and black lines indicate limits of 199 Monte Carlo simulation envelopes.
Spatial structure of Abies marocana 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 ag-gregation of saplings around females while mark analysis using k mm (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) identified mostly independence of tree dbh at greater distances using k mm (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 k mm (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 k mm (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 g 12 (r) or p 12 (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 higher Forest Systems August 2020 • Volume 29 • Issue 2 • e014 rate 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 indivi-duals, 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 fir. 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.