Tree mortality in mature temperate forests of central Mexico: a spatial approach

Tree mortality is an important process of forest stand dynamics and knowledge of it is fundamental to implement adequate management strategies. Subject to several factors, tree mortality can induce different spatial patterns on the remaining live and dead trees. While spatially clustered tree mortality in young forests is often driven by competition, in old-growth forests, spatially clustered tree mortality is often caused by disturbance agents. This study is focused on a spatiotemporal analysis of tree mortality in a mature temperate forest located in central Mexico dominated by Pinus montezumae and Alnus firmifolia. We used tree locations from a permanent plot (300 × 300 m) measured over a 20-year period. The results, from applying point pattern analysis, showed that the spatial pattern of all dead trees was clustered at short to medium distances, but showed no clear deviation from complete spatial randomness at longer distances. Similar results were found for P. montezumae and A. firmifolia. Using the bivariate mark-connection function (alive and dead trees), no tree mortality caused by competition was discernable, only A. firmifolia showed a tendency toward competition-introduced mortality around 15 m. Regarding forest structure, alive trees retained a clustered distribution and size heterogeneity at different distances during the measurement period. Thus, there was evidence that the resulting spatial pattern of tree mortality could be explained by disturbance agents such as droughts rather than tree competition. Therefore, the results of this study can contribute to implement management strategies based on the principles of continuous cover forestry and provide novel information regarding tree mortality in Mexican montane forests.


Introduction
Tree mortality is an essential process in forest ecosystems and knowledge of it is fundamental for sustainable forest management. It has multiple effects on forest structure Communicated by Lauri Mehtätalo.
1 3 and ecosystem functions. For instance, dead trees lead to open spaces that promote regeneration, incorporate nutrients into the soil and provide habitat for wildlife (Franklin et al. 1987;Landsberg and Sands 2011). Although many efforts have been made to comprehend its role on forest ecosystems, less attention has been focused on the spatial aspects of tree mortality (Larson et al. 2015). However, understanding the drivers of tree mortality and their effects on spatial patterns could be one of the key aspects for the adequate use, conservation, and restoration of forest ecosystems (Franklin et al. 1987;Stephens et al. 2008;Larson et al. 2015;Szmyt and Tarasiuk 2018).
In temperate forests, tree mortality is primarily influenced by competition, age, disease, and disturbance agents. Competition takes place especially in the early stages of forest development, when young trees intensively compete for sunlight, nutrients, and water (Peet and Christensen 1987; Landsberg and Sands 2011). Age-related tree mortality is usually determined by the forest succession and its evolution over time (Franklin et al. 1987). Finally, depending on the regime (Pickett and White, 1985), disturbances can dramatically induce tree mortality and reshape the spatial patterns in forest ecosystems; for example, forest fires, hurricanes, insect attacks, and severe droughts (Franklin et al. 1987;Agee 1993;Rodríguez-Trejo 2014;Choat et al. 2018).
Regardless of the ecosystem and causes, spatial patterns of tree mortality tend to be clustered (Kenkel 1988;Larson et al. 2015). However, reasons for this patterning differ between young and old-growth forests. Since mortality in young forests is more likely attributed to a competitivebased mortality such as self-thinning (Kenkel et al. 1997;Kenkel 1988;Wiegand and Moloney 2004), tree mortality in older forests is often related to disturbance agents such as wildfires (Franklin et al. 2002;Stephens et al. 2008). Therefore, although the spatial patterns of dead trees sometimes do not differ from young to old, a shift from competitive mortality in young to a non-competitive tree mortality in old-growth forests is expected (Franklin and Van Pelt 2004;Getzin et al. 2006;Larson et al. 2015).
Generally, the spatial relationship between living and dead trees can be described as spatial attraction or segregation. Attraction describes situations in which live and dead trees are spatially close together, which is usually attributed to a mortality by competition (Kenkel 1988;Kenkel et al. 1997;Wiegand and Moloney 2014). Spatial segregation could be present when dead and alive trees are spatially separated from each other. This might be explained by the agent-based tree mortality (Franklin et al. 1987;Landsberg and Sands 2011). In addition, spatial segregation is often related to the scramble competition, resulting in spatial patterns of clumped dead trees of the same species with more isolated surviving trees (Wiegand and Moloney 2014).
The spatial patterns of the remaining living trees are related to the processes described above. Thus, there is the hypothesis that the outcome of competitive-based mortality results in a more regular pattern of trees (Kenkel et al. 1997) and can be specially found in young forests (Larson et al. 2015). Contrastingly, in old-growth forests, agentbased mortality leads to spatial clustering in the remaining trees (Franklin and Van Pelt 2004). This clustering, together with size heterogeneity, has been connected to more resilient forests, especially to forest fires (Stephens et al. 2008). This provides important insights about the forest dynamics, that could be implemented into management strategies, grounded in the expectations of global change to come (Choat et al. 2018).
Like most ecological processes and spatial patterns, tree mortality is also likely to be scale-dependent. Point pattern analysis is a powerful tool for inferring underlying processes from the resulting scale-dependent patterns (Brown et al. 2011). Thereby, the main assumption is that the spatial patterns contain information about the processes that shaped them. By describing, analyzing, and modeling the spatial point patterns, i.e., the spatially explicit locations of all individuals in the study area, point pattern analysis allows exploring and confirming ecological hypotheses of the acting processes (Law et al. 2009;Brown et al. 2011;Wiegand and Moloney 2014).
The aim of this study was to deepen current knowledge of tree mortality in mature temperate forests. For this purpose, we used data from a permanent plot (9 hectares) located in central Mexico. The point patterns of all individual trees were repeatedly measured over 20 years. To the best of our knowledge, in this region, there is no study on spatial patterns of tree mortality, conducted at these scales, both spatially and temporally. Regarding the study area, our research questions were (i) Are dead trees completely spatially randomly distributed? (ii) Is tree mortality caused by competition? iii) Does tree mortality lead to spatial regular distribution? (iv) Does tree mortality lead to spatial size homogeneity? The answer to these questions could help to understand how this type of forest develops, and to utilize this information to implement the most appropriate methods for the use, conservation, and restoration of temperate forest ecosystems.

Study area
This study was carried out in the "San Juan Tetla" old research station (SJTRS), Chiautzingo, Puebla, Mexico (Fig. 1). The SJTRS used to be located besides the western part of the Ixta-Popo National Park, between the geographic coordinates of longitude [98.55°, 98.58°] and latitude [19.17°, 19.22°]. The soils are classified as Ando, or Humic Allophane with a mean slope between 30 and 40%, mainly facing the East exposition. The mean temperature is 9°C with a minimum of − 5°C and the precipitation is mainly concentrated from May to October, with 849.6 mm, which means 84.2% of the total, and from November to April 159.4 mm. The maximum monthly rainfall is in June, with 215.7 mm and the minimum in December, with 11.9 mm (Acosta-Mireles et al. 1997;Zepeda-Bautista and Acosta-Mireles 2000).
The vegetation type is classified as temperate forest with predominance of Pinus montezumae Lamb. This species is distributed in an elevation range from 2700 to 3450 m over the sea level, and it is commonly found in association with Pinus teocote Schltdl. & Cham., Pinus ayacahuite C. Ehrenb. ex Schltdl., Pinus leiophylla Schiede ex Schltdl. & Cham., Abies religiosa (Kunth) Schltdl. & Cham. and broad-leaved trees like Alnus firmifolia Fernald and Quercus sp. (Zepeda-Bautista and Acosta-Mireles 2000). This ecosystem has a potential fire regime of low severity frequent fires (Rodríguez-Trejo 2008;Jardel-Peláez et al. 2009), and P. montezumae is highly adapted to them (Rodríguez-Trejo and Fulé 2003). Nevertheless, it has been reported that fires in combination with diseases can increase the risk of tree mortality (Fonseca-González et al. 2016).

Data collection
The data were collected from a permanent plot (300 × 300 m, 9 ha in total) located in a mature forest, with three censuses over 20 years. The average tree-age was 70.5 years in the beginning of the measurement period (Zepeda-Bautista and Acosta-Mireles 2000). The plot was established and the first census was carried out by the National Institute of Forestry, Agriculture and Livestock Research (INIFAP, Spanish abbr.) in 1974. The second and third census were carried out in 1989-90 and 1995-96. The original plot consisted of 16 ha, however, we excluded some areas to maintain data consistency between 1974, 1989 and 1995. For more details about the data collection, see Zepeda-Bautista and Acosta-Mireles (2000). The censuses included all trees with a diameter at breast height (d) larger than 7.5 cm. All trees were tagged, the x and y coordinates recorded, and the following inventory variables were taken: tree diameter (d, cm), total height (h, m), crown width (cw, m), crown base height (cb, m). In addition, qualitative data, such as damages and condition (living or dead), as well as the species, were included.

Spatial distribution of dead trees: pair-correlation function g(r)
Spatial point pattern analysis uses the spatial location of all individuals within a study area described as points. To characterize spatial tree mortality patterns, we used several so-called second-order summary functions. These functions describe a point pattern based on its distances between points, allowing to describe the pattern at several scales or distances. We assumed that the underlying point process is stationary and isotropic, i.e., translation and rotation invariant. Thus, the mean number of points per area, the intensity lambda λ(x), is constant. All second-order summary statistics we restricted to a maximum distance between points of 18 m which equaled the assumed maximum interaction distance based on the largest crown diameter of 8.7 m. For all summary statistics used in this study, we applied the isotropic edge correction (Ripley 1977;Stoyan and Stoyan 1994).
First, we used the pair-correlation function g(r) (Stoyan and Stoyan 1994) which describes a point pattern by counting the expected density of points at a certain distance r of an arbitrary point of the point process. The g(r) function is defined as where K'(r) is the derivative of the Ripley's K-function (Ripley 1977). Thus, while Ripley's K-function is cumulative, i.e., using the expected density of points within a certain distance r, the pair-correlation function g(r) is non-cumulative, i.e., using the expected density of points at a certain distance r. This has the advantage that g(r) values at larger distances r are not influenced by smaller distance r and while g(r) = 1 (1) g(r) = K � r∕(2 r) indicates complete spatial randomness (CSR) at distance r, g(r) > 1 indicates that the points are clustered at distance r, whereas g(r) < 1 indicates regularity at distance r, respectively (Stoyan and Stoyan 1994;Wiegand and Moloney 2014).

Mark-connection function p ij (r)
Additionally to the spatial location, points can also have additional information related to them, so-called marks (Wiegand and Moloney 2014). Thus, we used the bivariate mark-connection function p ij (r) as summary function to describe the relationship between two discrete marks, namely dead and alive trees. For two marked points separated by a distance r, p ij (r) describes the probability that the first point has the mark i (i.e., dead) and the second point has the mark j (i.e., alive) (Illian et al. 2008;Wiegand and Moloney 2014). The mark-connection function p ij (r) is defined as In this case, the observed p ij (r) above (p ij (r) > p i p j ) and below (p ij (r) < p i p j ) the theoretical value indicates association or segregation, respectively (Stoyan and Stoyan 1994;Illian et al. 2008;Wiegand and Moloney 2014).

Mark-correlation function k mm (r)
Finally, to test the spatial correlation between live tree sizes, we used the mark-correlation function k mm (r). We used the tree diameter as size mark because generally the diameter is well correlated with tree height or tree crown, yet the diameter is easier to measure with less errors. The mark-correlation function describes the spatial correlation between continuous marks of two points (e.g., tree diameter) separated by distance r by comparing the mean of the marks in a pair of points to the mean mark based on all points (Stoyan and Stoyan 1994;Illian et al. 2008;Velázquez et al. 2016). The function k mm (r) can be estimated as where o is the origin point and r is a point at distance r, m(o) and m(r) are the marks assigned to them, and 2 is the squared mean mark as normalization factor. Thus, k mm (r) = 1 indicates an absence of correlation, values k mm (r) > 1 indicate a positive correlation and values k mm (r) < 1 indicate a negative correlation at a distance r (Stoyan and Stoyan 1994;Wälder and Wälder 2008).

Null models and simulation envelopes
In spatial point pattern analysis, null models act as in standard statistical null hypothesis testing and can be used in an exploratory or confirmatory purposes (Wiegand and Moloney 2014). In this context, null models usually use point process models to generate random patterns based on certain ecological hypotheses. Selection of an appropriate null model is one of the most challenging tasks (Wiegand and Moloney 2014) and depends on the specific question.
In this study, we used three null models. First, the simplest used null model was complete spatial randomness (CSR) for univariate patterns. Here, we used a homogeneous Poisson process to create random null model data. The homogeneous Poisson process is characterized by a constant intensity λ(x) throughout the study plot and no interaction between individual points. Second, we used a Thomas cluster process to simulate a clustered pattern. The Thomas cluster process simulates cluster centers using a Poisson process. Following, at the cluster centers (i.e., parent points), so-called offspring points are placed, whereas the number of daughter points is Poisson distributed and the locations are independent according to a two-dimensional isotropic Gaussian distribution (Illian et al. 2008). We fitted the Thomas cluster process to data that showed significant clustering only. Third, we used a random thinning null model to test if mortality leads into clustered point patterns. For this, we used the initial point pattern in the year 1974 and randomly removed points (i.e., simulated random mortality). The probability of each point to die corresponded to the mortality rates recorded over the measurement period of 1974 to 1995. Last, we used random labeling as null model. This null model differs from the previous used null models because it is not a point process model in a sense that point locations are simulated. Random labeling rather randomizes the marks of a bivariate pattern while keeping the point locations as they are (Wiegand and Moloney 2014). Random labeling is the appropriate null model if not the point locations are of interest, but a process acting on the already established locations, for example to test the random mortality hypothesis (a posteriori effects, Goreaud and Pélissier 2003).
To test the null hypothesis, we used Monte Carlo envelopes built from 199 replicates of the null models. From these simulations, the fifth largest and lowest values, respectively, were used to construct the simulation envelopes. For distances r at which the observed summary function is located within the simulation envelope, the null hypothesis is accepted (e.g., no significant departure from the CSR); if the observed summary function is located above the simulation envelope clustering/positive correlation is indicated, and correspondingly if the summary function is located below the simulation envelope regularity/negative correlation is indicated.
Because simulation envelopes test the null hypothesis at several distances r simultaneously, the test is susceptible to a Type I error inflation. Thus, simulation envelopes are no confidence intervals and do not allow to reject a null hypothesis with a certain significance level if the observed summary function is located outside the envelopes at a specific distance r. Nevertheless, they are still a useful tool because they show the potential result if a specific distance r would have been fixed previously (Baddeley 2014). To account for the multiple testing issue, we used a goodnessof-fit (GoF) test based on Loosmore and Ford (2006). The Diggle-Cressie-Loosmore-Ford test (DCLF test) uses the integral of the squared difference between the observed and simulated summary functions as reduced, single test statistic and consequently allows to accept or reject a null hypothesis (Loosmore andFord 2006, Baddeley 2014). The DCLF test was applied to the same distances as all secondorder summary functions, namely 0-18 m. Furthermore, because the observed summary functions could be located either above or below the simulated summary functions, we used a two-sided test. However, because the null models depend on parameters estimated from the data, the DCLF tests are conservative, i.e., the true significance level might be smaller than the reported significance level resulting in a wrong acceptance of the null hypothesis (Baddeley 2014). Thus, we additionally calculated all DCLF tests using a different summary function than used for the simulation envelopes, namely Besag's transformation of Ripley's K-function (Besag 1970). Nevertheless, no significant values changed, and we report all DCLF results based on the original secondorder summary functions. An alternative approach to the GoF test is global simulation envelopes which use a constant width based on the most extreme deviation from the theoretical null model value (Baddeley et al. 2014). All the analyses in this research were processed through R 3.6.0 (R Core Team 2021) using the spatstat R package (Baddeley et al. 2015). For specific questions, methods, and models, see Table 1.

Dead and alive trees tracking
Two species dominated the point pattern of dead and alive trees, namely P. montezumae and A. firmifolia. During 1974During -1989, the mean annual death was 17.5 trees, of these 7.7 were P. montezumae and 9.7 were A. firmifolia. However, by 1995, these numbers changed, mean annual death was 33 trees from 1989 to 1995, reducing the proportion of dead A. firmifolia, while it increased for P. montezumae (Table 2). It should be noted that this period was shorter (6 years), indicating that the death rate was higher in the case of P. montezumae.
Most of the trees that died, regardless of the species and period, had a rather small d (Fig. 2). However, while P. montezumae dead trees were found in all diameter categories in both measurements (1989 and 1995), for A. firmifolia, most of dead trees were found in the smaller categories. Besides, most of A. firmifolia larger dead trees were found by 1989 (Fig. 2).
The sizes of the alive trees by 1974 and 1995 showed a bimodal diameter distribution (Fig. 3). In both cases, the smaller diameters mostly belonged to A. firmifolia trees, while larger diameters mostly belonged to P. montezumae, meaning that the pines were the dominant trees in the study area. The shape of the diameter distributions had not changed remarkably in 15 years of measurements and only shifted from one diameter category to another, following a diameter increment given by the age (Fig. 3). Tree mortality leads to spatial regular distribution in young forests and clustered distribution in old-growth Spatial distribution of alive trees before and after the tree mortality period. Univariate pair-correlation function g(r) CSR null model Does the tree mortality lead to spatial size homogeneity of the remaining alive trees?
Spatially tree size heterogeneity is found in old-growth forests and related to natural perturbation regimes Spatial and size related distribution of alive trees before and after the tree mortality period. Mark-correlation function k mm (r) Random labeling null model

Spatial patterns of dead trees
The mapped data of all dead trees during the entire 20-year period are illustrated in Fig. 4. The overall pair-correlation function showed that tree mortality significantly differed from CSR at various distance after 20 years. All dead trees were aggregated at short to medium distances (0.0-18.0 m) (Fig. 5). The analysis by species showed some differences: The aggregation was significant at shorter distances for P. montezumae (< 12 m), while A. firmifolia was aggregated at almost all distances (Fig. 5). Simultaneously, the spatial point pattern of all alive trees in 1974 tended toward CSR for most spatial distances r and was at least less clustered for P. montezumae and A. firmifolia, respectively. Lastly, the random thinning null model resulted in less clustered point patterns for all trees as well as the point patterns separated by species (Appendix Figs. 1-7). Since the pair-correlation function g(r) revealed significant clustering of dead trees, we fitted the Thomas process to all dead trees. The Thomas process could model the observed pattern quite well and the pair-correlation function of the observed pattern was within the simulation envelopes at almost all distances r. The parameters of the final Thomas process model for the all dead trees group resulted in a density of 245.55 clusters in the whole area (26.8 ha −1 ), with a mean of 1.85 offspring points around parents points. P. montezumae showed a density of 199.16 clusters (21.74 ha −1 ) and a mean of 1.16 offspring points, while A. firmifolia showed 107.79 clusters (11.77 ha −1 ) with a mean of 2.07 offspring points (Fig. 6).

Mark-connection function p ij (r)
The testing for density dependencies between dead and alive trees p ij (r) showed no tendencies toward competition for most of the distances in the overall analysis (Fig. 7, left). However, this relationship showed segregation at short distances. The analysis by species indicated that the relationship for dead and all alive trees of P. montezumae had a tendency toward segregation around 12 m (Fig. 7, middle). However, the p value of the DCLF test was not significant. Nevertheless, for A. firmifolia, dead trees showed segregation from alive trees at short distances and a tendency toward attraction only significant around 15 m (Fig. 7, right).

Spatial patterns of alive trees: Pair-correlation function g(r)
In order to assess the changes in the spatial distribution of the remaining alive trees in 1995 in contrast to the initial alive trees in 1974, we made size classes based on percentiles (small = [d < = 33rd], medium = [> 33rd d < = 66th], and large = [d > 66th]). This provides consistency by using standardized tree diameter classes to each measurement period (1974 and 1995) and allows to have enough points to perform the spatial analysis. For the first measurement period (1974), the small class included trees with d < = 19 cm, the medium class trees with 19 < d < = 46.7 cm, and the larger class trees with d > 46.7 cm. For the last measurement period (1995), the small class included trees with d < = 28.7 cm, the medium class trees with 28.7 < d < = 56.6 cm, and the larger class trees with d > 56.6 cm. Spatial patterns for all alive trees were different among size classes, but remain almost the same through time. By 1974, the small alive trees were more aggregated than the other diameter classes in almost all distances, while the medium size class trees were aggregated at short distances and the larger trees were aggregated only at very short ones (Fig. 8, upper plots). After 20 years, in 1995, the remaining alive trees showed almost no difference in spatial distributions related to size classes. In a general analysis, the comparison of all the initial alive in 1974 vs. all the remaining alive trees showed that the spatial pattern slightly changes from clustered toward CSR (Appendix Fig. 8). Thus, this trend might continue in the future and if competition is strong enough even a regular pattern could occur (Couteron and Kokou 1997;Hesselbarth et al. 2018).
The Thomas cluster process was fitted only to the small and medium size classes for the alive trees both in 1974 and 1995. The parameters of the final Thomas process model for the small class [1974] group resulted in a density of 154.74 clusters in the whole area (17.03 ha −1 ), with a mean of 5.1 offspring points. The medium class [1974] showed a density of 940.89 clusters (102. 62 ha −1 ) and a mean of 0.84 offspring points. Finally, the small class [1995] has a density of 160.66 clusters (17.66 ha −1 ) with a mean of 3.94 offspring points and the medium class [1995] with a density of 727.41 clusters (79.52 ha −1 ) and a mean of cluster size  (Fig. 9). Thus, for the medium class [1974] and the medium class [1995], the Thomas cluster process was not an appropriate null model to describe to underlying point process because the mean number of offspring points was less than one point per cluster and the clustering was probably driven by other mechanism. However, also a heterogeneous null model did not result in a satisfactory fit.

Spatial patterns of alive trees: Mark-correlation function k mm (r)
The results of the Mark-correlation Function (k mm (r)) showed that the tree mortality did not affect the spatial size original distribution. In general, the analysis did not show a strong correlation between the marks; however, the alive trees at the beginning of the period had a tendency toward a negative correlation at different distances (6.0-18.5, 37.5-42.5 m), indicating spatial size heterogeneity, i.e., large and small trees being close together. This tendency

Fig. 6
Pair-correlation function g(r) of the fitted Thomas cluster process for the dead trees. The black solid line describes the observed g(r); the gray solid lines are the fifth largest and lowest values of the point-wise simulation envelopes based on 199 simulations of the fitted Thomas process; λ p and σ are the intensity and standard deviation estimated from the observed data was retained by 1995, only a few small differences at very specific distances were detected (Fig. 10).

Discussion
The focus of this research was to deepen current knowledge of tree mortality in mature temperate forests. According to our major findings, the spatial distribution of tree mortality i) occurred in clusters at different scales; ii) in general, did not show strong evidence of mortality by competition, iii) the spatial distribution of the alive trees did not show strong changes, and iv) the spatial size heterogeneity of the alive trees was maintained the entire analyzed period. These results are consistent with other studies carried out in mature and old-growth forests (Getzin et al. 2006;Larson et al. 2015).

Spatial patterns of dead trees
Trees that died in the study area and during the sampling term were spatially clustered at short and medium distances (0.0-18.0 m). In general, this clumped arrangement of dead trees is driven by mortality due to competition, often related to the self-thinning process in young forests (Kenkel 1988;Kenkel et al. 1997;Landsberg and Sands 2011;Wiegand and Moloney 2014) and by non-competitive or agent-based tree mortality, where mortality is caused by a variety of disturbance agents (White and Pickett 1985;Franklin et al. 2002;Landsberg and Sands 2011). In the study region, the most common disturbance agents are fire, diseases, and anthropogenic interventions (Rodríguez-Trejo 2008;Fonseca-González et al. 2016). However, by the time of this study, only a small fire was detected in the southern part of the plot that was contained before it caused a mayor damage to the trees.
In general, other causes of agent-based tree mortality are related to environmental factors, such as high temperatures and droughts (Franklin et al. 1987;Landsberg and Sands 2011;Wang et al. 2012). Throughout the region, there were severe drought seasons in 1983 and 1988 [− 3/− 4 values of the Palmer Drought Severity Index (PDSI)] (Stahle et al. 2016) before the first re-measurement in 1989, overlapping with El Niño/Southern Oscillation events (Wolter and Fig. 9 Thomas cluster process applied to the pair-correlation function g(r) of alive trees by class, where the black solid line describes the observed g(r); the gray solid lines are the 95% the fifth largest and lowest values of the point-wise simulation envelopes using 199 Monte Carlo randomizations of the Thomas cluster process; λ p and σ are the intensity and standard deviation from the fitted point process model Timlin 2011). Severe droughts cause stress on trees and can increase their probabilities to die (Franklin et al. 1987;Wang et al. 2012;Choat et al. 2018). This provides another perspective that might help explain the spatial patterns of dead trees in Fig. 5 and the mortality rates shown in Table 2 and Fig. 2, similarly to the study of van Mantgem and Stephenson (2007) in a temperate forest.
In this context, more research is needed in order to understand the effect of droughts on the species studied and to answer possible questions related to drought-induced tree mortality. For example, are these species tolerant or intolerant to drought and to what kind of drought are they adapted better. In this sense, McDowell et al. (2008) provide a theory on the hydraulic framework and the mechanisms of droughtinduced tree mortality. In their study, three mechanisms are hypothesized: (1) carbon starvation, referring to the reduction of photosynthesis caused by drought, (2) hydraulic failure, occurring when the drought induces desiccation and (3) biotic agents (insects and pathogens), referring to the synergy given by the combination of biotic agents and hydraulic failure or carbon starvation (McDowell et al. 2008).
In this respect, Himmelsbach et al. (2012) investigated the stomatal regulation of water status in a mixed pine-oak forest in northern Mexico. Where the species studied showed different water status regulation strategies, Pinus pseudostrobus was classified with isohydric water status regulation and Juniperus flaccida and Quercus canbyi with anisohydric water status regulation. This implies different responses to drought, e.g., while isohydric species respond better to severe but short drought events, anisohydric species are better adapted to long-term, low intensity drought events (Breshears et al. 2009;Himmelsbach et al. 2012). Defining these mechanisms is crucial to understand the drought-related tree mortality (McDowell et al. 2008) and the potential effects of climate change. Therefore, it is also important to establish long-term studies about plant strategies to drought such as the isohydric and anisohydric regulation (McDowell et al. 2008;Himmelsbach et al. 2012) and their effects on forest spatial structure.

Spatial relationship between dead and alive trees
Competitive tree mortality by suppression, i.e., alive trees surrounded by dead ones, was not supported by the data. Only A. firmifolia showed a tendency to aggregation with alive trees, although it only showed a deviation from CSR at distances of around 15 m. Therefore, it would not be appropriate to assume tree mortality related to competition, because the mean crown radius for A. firmifolia is about 2 m and for P. montezumae is about 3.5 m, meaning a low probability of crown overlapping. Other processes could be a possible explanation for this: i) age, as this species is a pioneer species (Gutiérrez-Guzmán et al. 2005) and whereas as the forest succession develops other species take its place, ii) disturbance agents, the severe droughts seasons in 1983 and 1988. Especially in the second case, severe droughts would exacerbate competition for water as a limited resource (Franklin et al. 1987;Landsberg and Sands 2011;Wang et al. 2012), resulting in clusters of dead trees, a situation that would be explained more by the scramble competition hypothesis (Wiegand and Moloney 2014).

Spatial patterns of alive trees
The uniform distribution theory, according to which tree mortality results in a regular spatial distribution of the remaining living trees, was not fully supported by the study. However, processes that usually lead to uniform distribution of the remaining trees are especially found in young forests (Kenkel et al. 1997;Larson et al. 2015). In other cases, oldgrowth forests showed no significant changes toward uniformity (Larson et al. 2015). Nevertheless, it is expected in old-growth forests that tree mortality by non-competitive agents leads to a heterogeneous distribution of alive trees (Franklin et al. 2002) and forests that have spatial and size heterogeneity become more resilient to disturbances (Stephens et al. 2008).
In this research, the rates of tree mortality maintained the size spatial heterogeneity through time (Fig. 10). This is to a certain extent, similar to fire-resilient forests of Jeffery pine in California (Stephens et al. 2008). This has important management implications, because historically in Mexico some common silvicultural treatments have been focused to produce uniform stands (Bray and Merino-Pérez 2005), contrasting with the natural structure of resilient forests. Besides, our results can be used as basis for designing forest restoration programs, based on spatial and size heterogeneity to promote stands close to natural structures similar to those found in this study.

Conclusions
To the best of our knowledge and especially for the Mexican forests, there is a lack of long-term data and spatially explicit tree mortality taken from medium-large plots. This study has both, a relatively large plot and standardized repeated measurements over 20 years; however, for this work, no spatial data of recruitment were taken. In addition, the minimal size of the trees to be included was 7.5 cm of diameter at 1.30 m above the ground. Therefore, we recommend that further research should include more information regarding seedlings and a standardized protocol for recording disturbances. Besides, although the measurement period was about 20 years, it is still a short period given the time needed for the forest to grow, thus we believe that an even longer time period is needed to have a better understanding of general forest processes.
Finally, this study provides relevant information on tree mortality and shows that it is a dynamic process that changes over time, space, and is species dependent. Understating tree mortality dynamics and the process that generate specific spatial patterns in forests is fundamental to the sustainable forest management (Franklin et al. 2002(Franklin et al. , 2007Stephens et al. 2008). In accordance with the principles of the continuous cover forestry or close to nature forestry (Gadow et al. 2002), it is necessary to understand how natural forests as dynamic systems develop through time and space in order to emulate these natural energy inputs and outputs throughout silvicultural treatments.
Author contributions EARC designed the study, conducted the data analysis, and wrote the manuscript; MHKH contributed to data analysis and writing; JGFG contributed to database pre-processing; MCM collected field information and generated the original database.
Funding Open Access funding enabled and organized by Projekt DEAL. The data collection was funded by the National Institute of Forestry, Agriculture and Livestock Research (INIFAP, spanish abbreviation). CONACYT and DAAD provided the first author with a joint grant (91680941) during his doctoral studies.

Data availability Not applicable.
Code availability The R code developed in this manuscript is available on request from the corresponding author.

Conflict of interest
We have no conflicts of interest to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.