Species habitat associations in an old-growth beech forest community organised by landslide disturbances

ABSTRACT In order to understand the processes that govern tree species’ spatial associations with habitats in landslide-affected forest communities, we investigated the habitat associations of six major tree species in an old-growth beech forest community located on a steep flank slope where landslides had occurred. All stems ≥ 5cm in diameter at breast height were mapped on a 1.14 ha polygonal plot and the topographic conditions (slope inclination and convexity), canopy state, and forest floor vegetation (i.e. dwarf bamboo and fern cover) were assessed. Most stems of Fagus crenata and Acer japonicum as well as many stems of Magnolia obovata belonging to the canopy layer were located on sites with low fern cover, whereas many stems of M. obovata below the canopy layer were associated with sites having high dwarf bamboo cover. The locations of the two Acer species belonging to the canopy and lowest layers coincided with sites having convex topography. Stems of Aesculus turbinata below the canopy layer were generally found on sites with gentle slopes, whereas most stems of Quercus crispula coincided with steeply sloping sites. While the results demonstrate that steep slopes created as landslide headscarps provide suitable habitats for species with less shade tolerance, it should be noted that there were few consistent patterns of habitat association across the layers. Our results suggest that unpredictable sporadic landslide disturbances with varying intensities and spatial scales, and their recursiveness, are at least partly responsible for habitat association patterns observed in the forest community.


Introduction
Natural disturbances play central roles in organizing forest communities (reviewed by Nakashizuka 2001). They impose spatial structure upon forest communities at various temporal and spatial scales, creating mosaics of vegetation patches at different regeneration stages (Thomas and Packham 2007). Spatial distribution patterns of trees after catastrophic disturbances, such as stand-initiating fires (Kenkel et al. 1997;Kashian et al. 2005), catastrophic typhoons or hurricanes (Yamashita et al. 2002;Xi et al. 2008) have been relatively well-studied. However, there is little published information on the spatial structure of well-mapped forest communities influenced by landslides (although see Martínez et al. 2010;Velázquez et al. 2014). Because many natural forests exist in mountainous regions with steep slopes such as those found across Japan, there is a need to understand how topographyrelated disturbances such as landslides affect the spatial structure of the forest communities in these regions.
A landslide is a mass movement on slopes; the term covers the both gradual slippage and more dramatic events that may occur at sites with steep gradients (Elias and Dias 2009;Namikawa et al. 2010). The original soil and vegetation is typically completely removed from the upper slopes of the disturbed area after a single landslide, creating areas with favourable light conditions that are amenable to colonisation by early successional or less shade-tolerant woody plants (Shiels et al. 2006;. Deposition zones on the lower parts of landslide slopes often have gentle gradients with varying micro-topographic relief and relatively humid soil conditions, which favours the establishment of tree species that require moist soil conditions (Lewis 1998;Koizumi 1999). Landslide disturbances also alter the forest floor vegetation of ferns, herbs, and bamboos ) that produce dense shade which prevents seedling establishment (Yamamoto et al. 1995;George and Bazzaz 1999). One would thus expect the spatial distributions of tree species to be structured by the formation of new habitats after landslide disturbances, generating species-specific habitat associations.
However, this simple model disregards several complexities. Landslides may occur repeatedly in certain areas with different timings, intensities, and spatial scales (Restrepo et al. 2009;. This will produce mosaic patterns in forest communities with individual patches being dominated by trees at different stages of regeneration (Seiwa et al. 2013;Velázquez et al. 2014).
Because not all landslides are catastrophic stand-initiating disturbances , they frequently cause changes in the forest floor environment (for instance by modifying the understory vegetation or soil moisture) without being so intense as to kill larger trees, i.e. trees belonging to upper structural layers [Okamoto and Higaki (2013), see also Fig. S1 of the Electronic Supplementary Material (ESM)]. Such disturbances will create a mismatch between the current forest floor environment and the conditions required by the dominant adult tree. Consequently, the spatial distributions of adult tree species are assumed to be inconsistent with their expected habitat associations. Therefore, to understand the processes that govern species' spatial associations with habitats in forest communities disturbed by multiple landslides, it is necessary to carefully analyse the effects of landslide disturbances on the spatial distributions of trees belonging to different structural layers.
Beech (Fagus crenata) forest is regarded as a climax forest community in the cool-temperate zone of Japan. There have been many studies on old-growth beech communities located on gentle slopes, demonstrating the effects of species-specific habitat associations on the organisation of tree communities (Yamamoto et al. 1995;Hirobe et al. 2015). There have also been studies on the species composition of beech forests subjected to landslide disturbances (Nakashizuka et al. 2003;Namikawa et al. 2010;Seiwa et al. 2013), and on successional patterns along the crosssection of the slope created by a single landslide (Koizumi 1999;Mishima et al. 2009). However, there have been little reported efforts to characterise the habitat associations of tree species in well-mapped beech forests that have been disturbed by multiple landslides differing in spatial and temporal scales.
The objective of this work is to bridge the knowledge gap discussed above by characterizing habitat associations of mapped trees of major species in an old-growth beech forest located on the lower flank of a steep flank slope that has experienced several landslides and exhibits considerable micro-topographical variation (Okamoto and Higaki 2013). The specific questions addressed in this work are: 1) Are there any tree species that exhibits a spatial association with habitats created by landslides? 2) Do the trends in species' spatial associations differ between structural layers in the studied area? In addition, we discuss ecological processes that could be responsible for the observed relationship between the spatial distributions of selected tree species and the landslide-affected habitats in the studied forest.

Study site
The studied stand was situated in an old-growth beech forest alongside the Sansuke-zawa creek, a tributary of the Ohkawa River in the Shirakami Mountains, located at the foot of Mt Yahougatake (40°30′ 7″ N, 140°10′ 48″ E; summit, 783 m above sea level [a.s.l.]). The Ohkawa River is one of the upper streams of the Iwaki River, which flows into the Sea of Japan. The forest soil and meteorological conditions in the Shirakami Mountains have been described in detail by Torimaru et al. (2013). In much of the stand, the understory was dominated by dwarf bamboo, Sasa kurilensis.
A field survey and geomorphological analysis indicated that the beech forests alongside the Sansuke-zawa creek were established on four major landslide blocks (Okamoto and Higaki 2013). This work focuses on a stand located on one of those landslide blocks, which is 120 m wide, 150 m long, and can be subdivided into several sub-blocks with widths ranging from 10 to 50 m (Fig. S2). The most recent landslide event in each block has been inferred to occur 30-180 years ago, implying heterogeneous disturbances at moderately large spatial scales (Okamoto and Higaki 2013). Thus, topographic variation in the study site has been influenced by past landslides.

Field methods
A 1.14 ha polygonal permanent plot was established in the studied stand in 2010. Because landslides are common in the Shirakami Mountains, the instability of the terrain often determines the area of easy accessibility, and this was the case with the polygonal plots considered here (e.g. Figure 1). All living woody stems with a diameter at breast height (dbh) of ≥ 5.0 cm were mapped to the nearest 0.1 m in three dimensions (i.e. xyz axes) using a TruPulse® 360 laser rangefinder (LaserTech Int., USA). We also recorded the species of each woody stem and measured its dbh to the nearest 0.1 cm. Using the categorisation proposed by Yamamoto et al. (1995), we identified three structural layers in the forest: layer Ithe canopy layer; layer IIthe subcanopy, comprising trees ≥ 8 m tall but below the canopy layer; and layer IIIthe understory layer, comprising trees < 8 m tall. Each stem was assigned to one of these layers based on its height.
In the autumn of 2010, prior to leaf-fall, canopy gaps were identified by visual inspection. To this end, we divided the entire plot into 5 m × 5 m quadrats, each of which was classified as having either a "closed canopy" (corresponding to canopy coverage of ≥ 30% at a height of >15 m) or a "canopy gap" (< 30% coverage). The cover of dwarf bamboo and that of ferns was estimated in each of 456 5 m × 5 m contiguous quadrats within the plot, and used to assign each quadrat to one of two cover-classes: high (≥ 30%) or low (< 30%). In addition to the three-dimensional coordinates of the base of every census tree, we determined the microrelief of the plot by measuring the elevation of 10 m × 10 m quadrats using the laser rangefinder.

Analysis of habitat association
We used the xyz coordinates of the census trees in 10 m × 10 m grid elements to estimate the elevation of 5 m × 5 m quadrats by linear interpolation. Then, to determine the local topographic state of each 5 m × 5 m quadrat, the index of convexity (IC) of Yamakura et al. (1995) was computed. The IC was defined as the difference between altitudes of the focal quadrat and the surrounding eight quadrats. Quadrats with positive and negative IC values have convex and concave slopes, respectively. The ground surface was approximated as a plane and the inclination of each quadrat was estimated based on the maximum inclination among those calculated for all pairs of corners of the focal quadrat.
To examine habitat associations of the tree species we used a generalised linear model (GLM), with a log link function and negative binomial error distribution to account for overdispersion in the response variable (number of stems of each species and layer in the 5 m × 5 m quadrats). The continuous values of convexity and inclination of the 5 m × 5 m quadrats were not directly measured in the field, but estimated by interpolation, so they may not be accurate enough for use as explanatory variables in habitat association analyses. Therefore, they were binarised to reduce unexpected artificially-introduced variations in the GLM analyses (see also Torimaru et al. 2013). Consequently, five environmental factors related to tree distribution, including two biotic factors (cover of dwarf bamboo or fern: high, 1; low, 0) and three abiotic factorsshading (closed canopy, 1; gap, 0), convexity (convex, 1; concave, 0), and slope inclination (steep, 1; gentle, 0), were used as predictor variables (see also Figure 1(a-e)). For each species and structural layer, Akaike Information Criterion (AIC) values were used to select the best candidate model. To evaluate the relative importance of each predictor, AIC values were computed for models based on all possible subsets of predictor variables. Using the procedure of Burnham and Anderson (2002), the difference in AIC value relative to the best candidate model (hereafter ΔAIC) was then calculated for each subset; any predictors included in all models with ΔAIC values ≤ 2 and < 4 were considered to have moderate and substantial effects, respectively, on the trees' heterogeneity of habitat association (see also Hirata et al. 2011). All GLM analyses were performed in R 3.1.2 (R development core team 2014). Because of the sample sizes (i.e. the number of trees) required for the habitat association analyses in each layer, we combined data from adjacent layers for A. turbinata (designated as layer I+II), and for M. obovata and Q. crispula (hereafter layer II+III) (see Results).

Floristic condition and physical properties in the stand
Twenty-seven tree and shrub species were observed in the 1.14 ha plot (Table S1 of ESM), which contained a total of 687 living stems with a total basal area of 38.9 m 2 ha −1 . This work focused on the six most abundant woody species at the study site, each of which individually accounted for at least 5% of the total density of living stems: Fagus crenata, Acer mono, Acer japonicum, Aesculus turbinata, Magnolia obovata and Quercus crispula. Those species collectively accounted for ca. 80.2% of the total stem density (Table 1, and see also Table S1 of the ESM).
The ground in the study plot slopes steeply from northwest to southeast, with its lowest point being in the southern corner of the plot. The maximum difference in elevation within the plot was 82.9 m. All 456 quadrats were classified as having either a convex (IC > 0 m, 210 quadrats) or concave (IC < 0 m, 246 quadrats) slope (Figure 1(a)). The mean gradient of the quadrats was 32.3°, with a maximum of 74.7°; their gradients ranged from gentle (6.0 -32.5°, 198 quadrats) to steep (32.5 -74.7°, 258 quadrats) (Figure 1(b)). Canopy gaps with varying sizes were identified over 0.13 ha (11%) of the plot (Figure 1(c)). Dwarf bamboo and ferns covered some of the forest floor: the proportion of 5m × 5m quadrats with coverage exceeding 30% was 31% and 51% for dwarf bamboo and ferns, respectively (Figure 1(d,e)). The evergreen fern Arachniodes standishii dominated in the ferncovered quadrats in the present study. This species was particularly abundant in the lower parts of the landslide deposits in the study region, where gravel and/or sand can be observed even at a depth of < 10cm (Mishima et al. 2009).

Associations with topographic types, canopy states, and/or forest floor vegetation
The spatial distributions of individual stems for each species (Figure 2) exhibited a relationship with the mapped topographic conditions and forest floor vegetation in some layers but not in the others (Figure 1 and Table 1). Negative coefficients of  Figure 2. Map of spatial distribution of individual stems of the six major tree species in the 1.14 ha plot. Large circles, medium-sized circles, and dots represent stems in layers I, II, and III, respectively; see Table 1 for the definition of the layers.
regression with fern cover were obtained for most stems of F. crenata (layers I and III) and A. japonicum (layer III), and approximately half of stems of M. obovata belonging to layer I ( Table 2), indicating that those stems were located on sites with low fern cover. Dwarf bamboo cover was identified as a predictor for the abundance of M. obovata stems in the lower layers (II+III), in which the stems were associated with sites that had high dwarf bamboo cover. Convexity was also a predictor for abundance of A. mono (layers I and III) and A. japonicum (layer III); positive coefficients indicating presence of their stems at sites with convex topography. The regression coefficients on steepness of slope were negative for A. turbinata (layer III), indicating an abundance of small adult trees on sites with gentle slopes, but positive for Q. crispula in both layers, indicating that most stems (all layers) occurred on steeply sloping sites. No significant associations with canopy state were found for any of the six major species.

Discussion
Spatial associations with habitats were detected in 11 out of 60 cases (Table 2). Only one species (Q. crispula) was associated with habitats created by landslides (steepness) across the layers, whereas no consistent patterns of habitat associations were found across layers for four species (F. crenata, A. mono, A. turbinata, and M. obovata). This contrasts with the result in an old-growth beech forests not subject to landslide disturbances in which the consistent patterns of habitat associations among layers were detected for four out of the five major species examined (Yamamoto et al. 1995). Slope inclination has previously been recognised as a key factor governing the spatial distribution of tree species (Ledo et al. 2013;Kitagawa et al. 2014;Sakai et al. 2014;Bourque and Bayat 2015). Most individuals of Q. crispula were found on sites with steep slopes that were evidently created by landslides (Fig. S3 of ESM); this is consistent with the conclusions of Yamamoto et al. (1995), who noted that this species can colonise areas of natural beech forests after landslide disturbances. Because Q. crispula is a gap-dependent species that is shade-intolerant or of intermediate shade tolerance (Masaki et al. 1992;Yoshida and Kamitani 2000;Kitao et al. 2006), the elimination of canopy trees during landslides could provide favourable light environments for its regeneration.
In contrast, small stems of A. turbinata were located at sites with gentle slopes, but no significant habitat associations were identified for this species in the upper layers (Table 2). In Japanese cool-temperate forests with landslides, Koizumi (1999) and Mishima et al. (2009) have found that adult trees of the species frequently occur on the lower parts of landslide deposits with gentle slope. The former study also found canopy trees on sites with steep slopes of scarps created by landslides, but the latter did not. Current landslide activities are very high in the site examined by Koizumi (1999) and our study site (Okamoto and Higaki 2013), but no recent landslides have been reported at the site examined by Mishima et al. (2009). This suggests that the site-specific landslide disturbance regimes affect the species' occurrences on sites with steep slopes, leading to variation in patterns of the species' habitat associations among layers.
Stems of M. obovata belonging to the lower layers were likely to be found in sites with dense Sasa cover, but this was not the case for those belonging to the canopy layer (Table 2), in contrast to patterns observed in previous studies of old-growth beech forests that are not subject to regular landslide disturbance (Yamamoto et al. 1995;Nakashizuka et al. 2003). The area considered here has been disturbed by several landslides with varying intensities and spatial scales (Okamoto and Higaki 2013), which may have eroded the Sasa cover over time, without killing the adult trees, in accordance with findings by Seiwa et al. (2013) that some adult trees could survive landslide disturbances.
The convex slopes in the studied plot could be classified into one of two types based on their contour lines (see Figure 1(a)): slopes of the main ridge extending from west to east in the center of the plot, and sites with microtopographic relief created by landslide deposits (on the eastern edge and the northern part of the plot, Figure 1 (a)). In A. mono, habitat association with convex slope was detected for the stems belonging to the lowest layers but not the upper layers (Table 2). Close examination of the spatial distribution map with contour lines shows that the species' stems were located not only on ridge slopes, Table 2. Regression coefficients of the generalised linear models with minimum AIC values for the six major woody species and layers in the 1.14 ha plot. Layers were defined in terms of crown position and tree height (see the footnote of Table 1 for details). Positive values indicate more frequent occurrence at the corresponding site, and negative the opposite. Bold and italic numbers denote variables that were included in all models with AIC differences (i.e. ΔAIC values) of < 4 and ≤ 2, respectively, compared to the best model. but also at sites with convex micro-topographic relief, especially at the eastern edge of the plot, irrespective of the layers (Figure 2). While the main ridge was created as a landslide block by an infrequent major landslide (see Study site), sites with convex micro-topographic relief might be created or eliminated by smaller, but more frequent landslides with varying intensities and spatial scales (Okamoto and Higaki 2013). These phenomena might be responsible for the variation in the patterns of habitat associations among layers for A. mono.
The sites with low fern cover were predominantly occupied by small stems of F. crenata and A. japonicum as well as canopy trees of F. crenata and M. obovata, but no consistent patterns of habitat associations with fern cover among layers were detected for the two canopy species in our study area (Table 2) (Koizumi 1999;Mishima et al. 2009). Generally, ferns comprise one of the most prominent plant life forms colonizing sites disturbed by landslides in temperate and tropical locations . However, these findings should be interpreted cautiously, as landslides have repeatedly occurred at some sites but less frequently, or not at all, in other parts of the study area (Okamoto and Higaki 2013). While one-off landslides allow ferns to persist after their colonisation, repeated landslides are likely to prevent their persistence. Thus, temporal scales of the emergence of habitats with low fern cover might highly vary, depending not only on the occurrence/absence of landslides but also on their recursiveness. This might partly explain the inconsistent patterns of habitat associations among layers.
In summary, in our landslide-influenced study area, landslides appear to alter the micro-topographic conditions and patterns of understory vegetation without killing the adult trees. In addition, the unpredictability of their intensity, spatial scale, timing, and recursiveness might contribute to the inconsistent patterns of habitat associations among structural layers of the tree species. Hence, we conclude that the landslides can, at least partly, cause deviation from expected species association patterns.