Supplementary Material to Souza 2017. Conifer Demography in Forest-Grassland Mosaics

Supplementary Material to Souza 2017. Conifer demography in forest-grassland mosaics: a landscape-scale study over a 24-year period. Botany.

Data collection

To develop a reliable identification practice of Araucaria canopies, a technician overlaid one of the two 2008 World View 1 images with stem maps of the species at 10 mapped 1-ha plots located at the São Francisco de Paula National Forest (29°25' S, 50°23' W). The National Forest is approximately 24 Km from the area covered by the aerial photographs. It was used as a training area for identification of Araucaria canopies since it harbours a long-term forest inventory and dynamics project with many mapped Araucaria trees. No further analyses were carried out with data from this area, which served exclusively for the training procedure described below. The National Forest is a 1606 ha reserve whose vegetation is composed of a mosaic of well-conserved patches alternated with forests that suffered selective logging of Araucaria and other timber species. Araucaria trees were mapped in the field in 2008 in ten 1-ha permanent plots, where all individuals with dbh ≥ 9.5 cm were measured for dbh and height. The plots correspond to a long-term project on forest structure and dynamics established in 2000 and monitored annually until 2009. Of the 10 plots, 5 were located in unlogged well-conserved forests, 3 were located in forests that suffered selective logging until 1945, and 2 were located in forests that suffered selective logging until 1989. In 2000 all trees in the 10 plots were mapped using a GPS receiver. Since the accuracy range of the GPS under the forest canopy was not great (approximately 3-5 m), the coordinates obtained this way were then checked in the field using measurements between each tree and edges of the corresponding plot (see Souza et al. (2008) for a detailed description of the study plots). The technician identified all Araucaria crowns in the World View image and checked whether these matched the adult Araucaria that were mapped in the field in the 10 ha. After some trial and error this matching was achieved, and the same technician then marked all Araucaria crowns in the 1984 aerial photographs and the 2008 World View image corresponding exactly to the aerial photographs.

Data analysis

Explanatory variables (distance from forest patch edge, patch area, number of first- and second-order conspecific neighbors) were subjected to exploratory data analyses (Zuur et al. 2010) in order to check for outliers, heterogeneity of variance, normality, and collinearity. Forest patch area was log-transformed in order to reduce a few extreme values. Variance inflation factors (VIFs; R function corvif from Zuur et al. 2010) calculated using all covariates indicated that collinearity was not an issue (all VIFs < 3.1). No significant interactions were found between the explanatory variables.

In order to choose the modeling approach that best accounted for the spatial structure in the data, I compared three different autocorrelation descriptors. First, I fitted a Generalized Linear Mixed Model (Zuur et al. 2009) to the data using the environmental variables as fixed factors and the affiliation of each individual Araucaria to either one of the two 1984 aerial photo mosaics as a random factor using the glmer function of the lme4 package (Bates et al. 2015). Photomosaic affiliation was explicitly included in the model in order to control for local history, soil type, biotic environment and other unknown site-specific differences (Zuur et al. 2009). Second, I fitted the data to autologistic or autonormal Generalized Linear Models (Besag 1974), which consist of including an additional term in the model (the autocovariate) to represent the influence of neighbouring observations, using the autocov_dist function of the spedp package. I used weights by the squared inverse distances of all data points in the neighborhood as the weighting scheme in the autocov_dist function (type = inverse.squared in the autocov_dist function), as suggested by Bardos (2015). Finally, I fitted a GLM to which I added up to fifth-degree polynomials of the coordinates as explanatory variables(Borcard et al. 2011), using R's base function glm. Note that the addition of polynomials to a GLM or a linear regression model do not make these models non-linear, it only accounts for spatial structure in the residuals. Spatial autocorrelation analysis of residuals from each model was performed using Moran's I correlograms as described previously, and the modeling approach that best controlled for spatial dependency was chosen. Here, I only report results from polynomial GLMs, since only these produced residuals with no spatial autocorrelation. In all cases, a full model was fitted first, followed by manual variable removal based on Akaike’s information criterion (AIC, Akaike 1974). Explanatory variables, polynomials included, were centered before analysis (Borcard et al. 2011).

Because my data consisted of maps of points, I used edge-corrected spatial point pattern analyses based on Ripley’s K-function (Ripley 1977), and its bivariate and marked extensions to test my sixth hypothesis as well as to assess the spatial structure of survivorship, crown growth and crown breakage and their associations with landscape features. These methods consist of density functions using the second moment (i.e., all point-to-point distances for a statistical description of two-dimensional distribution patterns). The calculation of K for a range of distances allows for the study of spatial pattern at a range of scales. Bivariate analyses of the spatial interaction (association/repulsion) between two groups of points (e.g. two species) were based on the K12 test, where distances are computed between points of the two groups (1 and 2), instead of points within the same group. Marked analyses are used when a numerical datum like size or growth is attached to each point in the pattern (Diggle et al. 1983). I used bivariate analyses to test for spatial association between trees in the grassland and trees in forest patches, and marked analyses to test for spatial structure in crown growth. Regarding bivariate analyses, as I studied the spatial interactions between two a priori different populations, the relevant null hypothesis is the one of independence and not of random labeling (Goreaud and Pélissier 2003). Univariate analyses tested the significance of the departure of observed samples from randomness. All spatial analyses were carried out with the ads package in R(Pélissier and Goreaud 2015).

My spatial data were divided in the two separated photomosaic areas and thus could not be readily analyzed. In order to convert the two sets of points into a single spatial object, we delimited a complex sampling window for spatial analyses consisting of two polygons (function swin of the ads package) corresponding to the shapes of the two aerial photo mosaics. This was accomplished by removing triangular surfaces from an initial sampling window of rectangular shape (function swin with option triangles set to an object containing triangles for exclusion, Fig. S1). The data of the two photo mosaics could thus be analyzed as a single point pattern. In point pattern analysis, edge effects arise from the fact that points near the edge of the study area tend to have higher nearest-neighbor distances, even though they might have neighbors outside of the study area that are closer than any inside it. Edge effects were corrected using Ripley’s unbiased local correction of edge effects for complex sampling windows (Pélissier and Goreaud 2015). Univariate, bivariate, and marked analyses were implemented with the kfun, k12fun, and kmfun functions, respectively. The derived variable L(r) (or L12(r) and gm(r) for the inter-type and marked analyses, respectively) enables the interpretation of the type of spatial pattern as a function of distance by plotting L(r) against r. For a completely random pattern L(r) = 0; L(r) becomes negative when the pattern is regular and positive when trees are clustered. Similarly, values of L12(r) = 0 indicate independence of the two types of points; L12(r) < 0 indicate repulsion effects, while L12(r) > 0 indicate attraction effects (Goreaud and Pélissier 2003).

A 10.0 m step was used for all point pattern analyses, which were computed up to 1000 m, a distance corresponding to half the smaller side of each aerial photo mosaic. Significance of departure from a random spatial pattern of null association between pairs of points was estimated by a Monte Carlo procedure to simulate randomly generated plots of the same density and dimensions as the observed plot. I used 1000 simulations of each analysis. In the bivariate analyses, new random coordinates were assigned to grassland trees only (the ‘dependent’ group). In marked point analyses, local Monte Carlo confidence limits were estimated at each distance r, after reallocating at random the values of the mark variable over all points of the pattern, the location of trees being kept unchanged (Pélissier and Goreaud 2015).

References

Akaike, H. 1974. A new look at the statistical model identification. IEEE Trans. Autom. Control 19: 716–723.

Bardos, D.C., Guillera-Arroita, G., and Wintle, B.A. 2015. Valid auto-models for spatially autocorrelated occupancy and abundance data. Methods Ecol. Evol. 6(10): 1137–1149. doi:10.1111/2041-210X.12402.

Bates, D., Machler, M., Bolker, B., and Walker, S. 2015. Fitting Linear Mixed-Effects Models Using lme4. J. Stat. Softw. 67(1): 1–48.

Besag, J. 1974. Spatial interaction and the statistical analysis of lattice systems. J. R. Stat. Soc. Ser. B 36: 192–236.

Borcard, D., Gillet, F., and Legendre, P. 2011. Numerical ecology with R. Springer, New York.

Diggle, P.J., Sibson, R., and Cohen, J.E. 1983. Statistical analysis of spatial point patterns. Academic Press, Lodon. Available from 243.

Goreaud, F., and Pélissier, R. 2003. Avoiding misinterpretation of biotic interactions with the intertype K 12 -function: population independence vs. random labelling hypotheses. J. Veg. Sci. 14(5): 681–692. doi:10.1111/j.1654-1103.2003.tb02200.x.

Pélissier, R., and Goreaud, F. 2015. ads Package for R: a fast unbiased implementation of the K-function family for studying spatial point patterns in irregular-shaped sampling windows. J. Stat. Softw. 63(6): 24.

Ripley, B. 1977. Modelling spatial patterns. J. R. Stat. Soc. B 39: 172–212.

Souza, A.F., Forgiarini, C., Longhi, S.J., and Brena, D.A. 2008. Regeneration patterns of a long-lived dominant conifer and the effects of logging in southern South America. Acta Oecologica 34: 221–232.

Zuur, A.F., Ieno, E.N., and Elphick, C.S. 2010. A protocol for data exploration to avoid common statistical problems. Methods Ecol. Evol. 1(1): 3–14. Available from

Zuur, A.F., Ieno, E.N., Walker, N., Saveliev, A.A., and Smith, G.M. 2009. Mixed effects models and extensions in ecology with R. Springer New York, New York, NY. doi:10.1007/978-0-387-87458-6.

Fig. S1. Spatial pattern of all mapped trees. Numbered triangles correspond to areas excluded from analyses with function swin (option triangles) of package ads.

Fig. S2. Spatial pattern of surviving and dying trees in the grassland and in forest patches.

Fig. S3. Spatial pattern of trees with crown breakage in the grassland and in forest patches. Circle sizes are proportional to individual Araucaria tree crown growth.

1

Table S1. Generalized Linear Models for selected predictor variables and survivorship of Araucaria angustifolia trees in forest patches. Models in bold are the best models.
Model / Fixed Variables / Spatial Variables / AIC / Spatial Autocorrelation of Residuals
Survivorship in forest patches
M1 / ICD+NCN1+NCN2+DFE+PA / None / 321.1 / Significant
M2 / DFE+PA / None / 319.9 / Significant
M3 / DFE+PA / X+X2+X3+Y+Y2+Y3+XY+X2Y+XY2 / 313.8 / Non-significant
M4 / DFE+PA / Y3+X2Y / 303.8 / Non-significant
Survivorship in the grassland
M5 / ICD+NCN1+NCN2 / None / 136.7 / Non-significant
M6 / NCN1 / None / 132.5 / Non-significant
Crown Growth in forest patches
M7 / ICD+NCN1+NCN2+DFE+PA / None / -8922.33 / Significant
M8 / ICD+NCN1+NCN2+DFE+PA / Full sixth-degree polynomials / -9086.78 / Significant
M9 / ICD / X2+X3+X5+Y+XY+Y2+XY3+Y5 / -9094.83 / Significant
Crown Shrinkage in the grassland
M10 / ICD+NCN1+NCN2 / None / 247.47 / Significant
M11 / ICD / X+X2+X3+Y+Y2+Y3+XY+X2Y+XY2 / 221.00 / Non-Significant
M12 / ICD / X+X2+Y2+Y3+XY+X2Y+XY2 / 217.39 / Non-significant

Continue.

Table S1. Continuation.

Model / Fixed Variables / Spatial Variables / AIC / Spatial Autocorrelation of Residuals
Crown Shrinkage in forest patches
M13 / ICD+NCN1+NCN2+DFE+PA / None / 999.78 / Significant
M14 / ICD+NCN1+NCN2+DFE+PA / Full fifth-degree polynomials / 895.46 / Non significant
M15 / ICD+NCN1+NCN2+PA / Full fifth-degree polynomials / 886.51 / Non significant
ICS = Initial Crown Diameter
NCN1 = First Order Number of Conspecific Neighbors
NCN2 = Second Order Number of Conspecific Neighbors
DFE = Distance from Forest Edge
PA = Patch Area

1

Table S2. Best Generalized Linear Model for selected predictor variables and survivorship, growth, and occurrence of crown shrinkage of Araucaria angustifolia trees.
Variable / Estimate / SE
Survivorship in forest patches
Intercept / 3.84 / 0.24
Distance from Forest Edge / 0.02 / 0.01
Patch Area (log) / 0.36 / 0.11
X3 / -3.38E-11 / 7.33E-12
X2Y / 2.35E-11 / 1.08E-11
Survivorship in the grassland
Intercep / 2.34 / 0.24
Number of First-Order Conspecific Neighbors / -0.54 / 0.24
Crown Growth in Forest Patches
Intercept / 0.0015 / 3.02E-05
Initial Crown Diameter / -0.0002 / 1.58E-05
X2 / -4.96E-07 / 1.41E-07
X3 / 4.65E-13 / 1.32E-13
X5 / -1.63E-25 / 4.65E-26
Y / 6.80E-01 / 2.80E-01
XY / 4.40E-08 / 1.25E-08
Y2 / -6.78E-08 / 2.74E-08
XY3 / -3.21E-22 / 9.12E-23
Y5 / 2.31E-29 / 8.74E-30
Crown Shrinkage in the grassland
Intercept / -0.21 / 0.20
Initial Crown Diameter / 0.47 / 0.09
X / 2041.00 / 841.60
X2 / -1.62E-04 / 6.51E-05
Y2 / 2.32E-05 / 9.83E-06
Y3 / -2.23E-12 / 9.67E-13
XY / -5.79E-04 / 2.47E-04
X2Y / 2.40E-11 / 9.63E-12
XY2 / 4.09E-11 / 1.82E-11

Continue.

Table S2. Continuation.

Variable / Estimate / SE
Crown Shrinkage in forest patches
Intercept / -1.76 / 0.12
Initial Crown Diameter / 0.70 / 0.05
Number of First Order Neighbors / 0.19 / 0.20
Number of Second Order Neighbors / -0.16 / 0.08
Patch Area (log) / 0.0014 / 0.0005
X / 1.09E+09 / 4.32E+08
X2 / -1009.00 / 408.40
X3 / -0.00002 / 0.00005
X4 / -1.03E-12 / 4.23E-12
X5 / -1.29E-19 / 1.93E-19
Y / 1.62E+08 / 5.48E+07
XY / -488.70 / 191.30
X2Y / 0.0005 / 0.0002
X3Y / 6.99E-12 / 1.60E-11
X4Y / 2.03E-19 / 6.54E-19
Y2 / -25.57 / 8.32
XY2 / 0.00007 / 0.00003
X2Y2 / -6.79E-11 / 2.65E-11
X3Y2 / -5.50E-19 / 1.20E-18
Y3 / 1.42E-06 / 4.35E-07
XY3 / -3.61E-12 / 1.40E-12
X2Y3 / 3.39E-18 / 1.31E-18
Y5 / -1.00E-21 / 2.78E-22

1