The elevational distribution of Fagus crenata regeneration at low elevations in the Shirakami Mountains, a World Natural Heritage site in Japan

ABSTRACT The Shirakami Mountains are a World Natural Heritage site because of their vast beech (Fagus crenata) forests and high diversity of flora and fauna. As global warming progresses, the area suitable for F. crenata growth is expected to shift from low to higher elevations. However, no studies have examined the current elevational distribution of the species’ regeneration. Here, we investigated the amount of F. crenata regeneration in the Shirakami Mountains and compared this amount among different elevations. Our nine study sites covered forests ranging from 269 m to 662 m, which covers the low to middle parts of the elevational distribution of F. crenata. At each site, we established three 20 m × 20 m plots and conducted a tree census. We approximated the regeneration density as a function of the diameter at breast height (DBH) distribution of F. crenata at each site using a power function: density = a DBHb. We found that a and b decreased and increased, respectively, as elevation decreased. Generalized linear mixed model analysis also indicated that decreasing elevation significantly decreased the number of seedlings of F. crenata. These results suggest decreasing regeneration at low elevation. Although it is difficult to identify the main causes of this pattern, we predict that regeneration may become even poorer if climate warming continues. Long-term monitoring along elevational clines will help us track the detailed effects of global warming on the beech forests.


Introduction
The Shirakami Mountains were designated a World Natural Heritage site in 1993 under criterion ix, "ecosystem," because of their vast forests of beech (Fagus crenata Blume) and the high diversity of flora and fauna.The forest covers 16 971 ha of Japan's Aomori and Akita prefectures and provides indispensable ecosystem services for local people, such as supplying many natural products, recharging water resources, and preventing surface erosion.Many animals live in the Shirakami Mountains, including notable species such as the serow (Capricornis crispus Temminck), Japanese dormouse (Glirulus japonicus Schinz), black woodpecker (Dryocopus martius Linnaeus), and golden eagle (Aquila chrysaetos Linnaeus).
Recent data suggest that the mean temperature in the Shirakami Mountains has been continuously increasing, accompanied by decreasing snowfall (Fig. S1).Because of this change, the area that is unsuitable for the growth of F. crenata is likely to be expanding at lower elevations.Indeed, many studies have shown that tree species are changing their distribution toward the poles and higher elevations (Walther et al. 2002;Penuelas et al. 2007;Beckage et al. 2008;Harsch et al. 2009;Feeley et al. 2011;Shimazaki et al. 2011;Koide et al. 2017Koide et al. , 2022;;Gatti et al. 2019).Wason and Dovciak (2017) suggested that distribution of American beech (Fagus grandifolia Ehrh.) is likely expanding upslope, although other factors, such as land-use changes, may mask species range shifts caused by changing climate.In Germany, the growth of European beech (Fagus sylvatica L.) has decreased at lower elevations but increased at higher elevations (Dulamsuren et al. 2017).In Japan, the relative abundance of evergreen broadleaved trees increased near the colder boundaries of their range, coinciding with a decrease in the relative abundance of deciduous broadleaved trees (Suzuki et al. 2015).A similar change was observed at the colder edge of the distribution of deciduous broadleaved forest communities, in which the abundance of deciduous trees is increasing, accompanied by decreasing abundance of boreal coniferous species.
According to the models of the potential impacts of climate change on canopy tree species in Japan, Matsui et al. (2018) predicted that the beech-dominated forest type will be replaced mostly by oak-dominated forest types.Their size-based species distribution model also implies that the potential habitat of juvenile F. crenata would be lost at warmer and snowless sites (Koide et al. 2016).In the Shirakami Mountains, a previous study based on ecological niche modeling predicted that the area suitable for F. crenata will decrease to nearly 0% within the World Natural Heritage area by 2100 as current global warming progresses (Matsui et al. 2007).These predictions suggest that F. crenata forest will gradually decline at the lower boundary of its distribution but will expand its range toward higher elevations.However, no studies in the Shirakami Mountains have examined the current elevational distribution of its regeneration.
Here, we investigated the amount of F. crenata regeneration in the Shirakami Mountains and compared this amount among different elevations.We paid particular attention to the lower boundary of F. crenata's distribution because that area may no longer be a suitable habitat for F. crenata in the near future due to current warming.In addition, we compared the pattern of regeneration with those of Magnolia obovata Thunb.and Acer pictum Thunb.subsp.mayrii, which are canopy tree species frequently associated with F. crenata in the Shirakami Mountains.The comparison allows us to understand whether the pattern observed in F. crenata is common to the other species.

Study species
Fagus crenata is one of the dominant tree species in the cool-temperate zone of the Japanese archipelago.It is a shade-tolerant, late-successional tree species and develops a significant seedling bank on the forest floor (Nakashizuka 1987;Koike 1988).Seedlings cannot grow rapidly under a closed canopy (Nakashizuka 1987) and have low survival rates under dwarf bamboo (Abe et al. 2002).Beech forest in Japan is unique due to its high understory coverage by dwarf bamboo.Thus, gap formation and the death of dwarf bamboo are important to promote F. crenata regeneration because these changes improve light conditions and decrease seed predation by rodents (Abe et al. 2001(Abe et al. , 2005)).Fagus crenata exhibits synchronized seed production, and its masting cycles affect the number of seed predators (Kon et al. 2005).

Study sites
The study was conducted in the northern part of the Shirakami Mountains in Aomori Prefecture, Japan (Figure 1).Beech forests in Japan are classified into the Japan Sea type and the Pacific Ocean type, which are distinguished primarily by maximum snow depth, with the separation at 50 cm (Shimano 2006).The Shirakami Mountains belong to the Japan Sea type, which has greater snowfall.
We selected nine study sites in the Shirakami Mountains and established three 20 m × 20 m plots at each site, for a total of 27 plots (Figure 1; Table S1).The elevation of the study sites ranges from 269 m asl (Lake Juniko site) to 662 m asl (Tsugaru Pass site).All plots were located in forests with closed canopies within the cool-temperate deciduous broadleaved forest zone.Matsui et al. (2007) found that the lower elevation boundary for beech forests in the Shirakami Mountains was 260 m, and the upper was 1070 m.Therefore, the elevation range of our study sites covers the low to middle parts of the elevational distribution of F. crenata in the Shirakami Mountains.Matsui et al. (2007) predicted the lower boundary of the potential distribution of F. crenata as 43 m based on the warmth index (Kira 1948).The area near this boundary has been disturbed by deforestation.Microtopography of the study sites was similar, with flat to gently sloped terrain.

Tree census
Our tree census was conducted from August to September 2021 at the 27 plots.We recorded the diameter at breast height (DBH, the diameter at a height of 1.3 m) of all trees in each plot.For individuals with DBH < 2 cm and ≥1.3 m in height, we recorded only the number of individuals.Hereafter, we call these individuals "saplings."In addition, we selected one corner of each 20 × 20 m plot described above and counted the number of seedlings of F. crenata by walking within the plot, along the edge, until the total area examined reached 25 m 2 (typically, a transect that covered 1 m × 25 m).Seedlings were defined as individuals shorter than 1.3 m in height.We determined the coverage by dwarf bamboo within each plot at 10% intervals since this species is known to affect F. crenata regeneration (Nakashizuka and Numata 1982).

Data analyses
We summarized the data from the tree census (see Tables S2a-i) and analyzed the relationships between the amount of F. crenata regeneration and elevation.First, we calculated the juvenile-to-canopy tree (J/C) ratio, where J represents the number of juvenile trees in the plot, which we defined as trees with DBH < 10 cm and ≥1.3 m in height, and C represents the number of canopy trees, which we defined as trees with DBH ≥30 cm (Shimano and Okitsu 1994).The ratio therefore indicates the relative proportion of juvenile trees of F. crenata.In this method, a higher J/C ratio means better regeneration.For this analysis, we pooled the data from the three plots at each site.The ratio was also calculated for M. obovata and A. pictum subsp.mayrii.
Second, we examined the shape of the DBH class distribution (Shimano 2000) of F. crenata.We constructed this distribution for each site by compiling the data from the three study plots with 5 cm class intervals.We transformed the data into a double logarithmic graph and calculated the coefficients a and b for the power function density = a DBH b .We subsequently compared the parameter values among the sites.When the proportion of small trees increases, a and b increase and decrease, respectively.In constructing the DBH class distribution, we included the sapling data in the lowest DBH class (i.e. 1 to 6 cm).We then quantified the relationships between the obtained values of a and b and the site's elevation by using a linear model.We did not conduct this analysis for M. obovata and A. pictum subsp.mayrii due to their small sample sizes.
Finally, we analyzed the relationships between elevation and the number of F. crenata seedlings and saplings by using a generalized linear mixed model (GLMM).In these analyses, the response variable was the number of seedlings or saplings, and the explanatory variables were elevation, coverage by dwarf bamboo, and the BA of F. crenata, excluding individuals with DBH <30 cm in each plot.BA of F. crenata with DBH ≥30 cm was included in the model as an indicator of the size of the seed source.We modeled the error distribution as a Poisson distribution with a log-link function.Site was treated as a random effect.The explanatory variables were standardized before the analyses.We selected the bestfit models using Akaike's information criterion (AIC); the model with the smallest AIC was selected as the best model.We also selected models with a similarly good fit (i.e.models whose AIC differed by <2 from the best model).For references, we performed the same analyses for M. obovata and A. pictum subsp.mayrii.For these species, we used the number of saplings in each plot as the response variable.We conducted the GLMM analyses using version 4.1.2 of the R software (R Core Team 2021) and the packages lme4 (https://cran.r-project.org/web/packages/lme4/index.html),lmerTest (https://cran.r-project.org/web/packages/lmerTest/index.html),and MuMIn (https://cran.r-project.org/web/packages/MuMIn/index.html).

Stand structure
In the tree census, F. crenata was the dominant tree species at all of our study sites; its BA was the largest at 22 out of 27 plots (Table S2).The average F. crenata BA was 31.90 m 2 /ha, and ranged from 13.67 m 2 /ha (TGB site, average of three plots) to 55.46 m 2 /ha (LJ site).Acer pictum subsp.mayrii and M. obovata had the second-and third-largest BAs, respectively.Other major associated tree species included Aesculus turbinata Blume, Quercus crispula Blume, and Tilia japonica (Miq.)Simonk.(Figure 2).At two of the lowelevation sites (LJ and BGSNT), the proportions of A. pictum subsp.mayrii and M. obovata were relatively high in the sapling population.The coverage by dwarf bamboo averaged 22%, with values ranging from 10% to 80% (Fig. S2).

Relationship between elevation and F. crenata regeneration
The F. crenata J/C ratio was relatively low at the low-elevation sites, although the relationship was not statistically significant (Figure 3a; p = 0.21).The highest J/C ratio (15.7) was recorded at the TGB site (elevation 559 m).In our analysis of the shape of the DBH class distribution for F. crenata, the a parameter increased slowly with increasing elevation, but the increase was only marginally significant (R 2 = 0.42; p = 0.06; Figure 4, Fig. S3).In contrast, the b parameter decreased significantly with increasing elevation (R 2 = 0.47; p = 0.04).
According to the GLMM analyses, the number of F. crenata seedlings was best explained by the model that included elevation and the BA of F. crenata canopy trees (Table 1).The full model that included elevation, coverage by dwarf bamboo, and the BA of F. crenata canopy trees provided a comparably good fit (i.e.AIC differed by <2 units).For the number of F. crenata saplings, the two models with comparable values of AIC provided comparably good fits: both included coverage by dwarf bamboo (significantly negative) and the BA of F. crenata canopy trees (significantly negative) as explanatory variables.Elevation was also included in the second model.

Comparisons of the regeneration patterns with M. obovata and A. pictum subsp. mayrii
The J/C ratio of M. obovata (Figure 3b) showed a similar trend to that of F. crenata, but there was a marginally significant positive relationship (p = 0.10).That is, the J/C ratio decreased with decreasing elevation.In contrast, A. pictum subsp.mayrii showed no significant relationship with elevation (Figure 3c; p = 0.44).In the GLMM analysis, elevation was not significantly related to sapling numbers for M. obovata and A. pictum subsp.mayrii (Table 1).However, the number of M. obovata saplings was significantly positively related to the BA of canopy trees, and the number of A. pictum subsp.mayrii saplings was significantly negatively related to the coverage by dwarf bamboo and to the BA of canopy trees.

Discussion
In the present study, we examined the relationship between elevation and regeneration of F. crenata.In the GLMM analysis, elevation was selected in one or several best-fit models as an explanatory variable for the number of seedlings and saplings of F. crenata (Table 1).Fagus crenata regeneration was limited at lower elevations as compared to that at middle elevation range, suggesting a clear effect of elevation on the number of seedlings.Moreover, the power analysis of DBH-class distribution indicates that parameters a and b decrease and increase, respectively, with decreasing elevation (Figure 4), though the increasing trend of parameter a was marginally significant.Parameters a and b represent tree density and decreasing density when tree size is doubled, respectively.Thus, the results suggest that the tree density was higher and the density increased more sharply as DBH class decreased at high elevation sites than that at low elevation.Shimano (2006) demonstrated that the population parameters of F. crenata differed markedly between the Pacific and Japan Sea populations; the latter had much better regeneration than the former.Our data demonstrated that elevational differences also existed in the forests of the Shirakami Mountains.Additionally, the J/C ratio tended to be low at the low-elevation sites (Figure 3) although the relationship was not statistically significant due to the high variation.
The elevational pattern of plant distribution and associated changes of species composition have been studied for tree species in other countries.In a Spanish study (Penuelas et al. 2007), the relative abundance of European beech (F.sylvatica) at the lower limit of the vertical distribution was half that at the center and upper limit of the species' distribution.The authors suggested that the result was likely due to current global warming.Furthermore, holm oak (Quercus ilex L.) in small size classes showed very high recruitment at the lower limit of the distribution of F. sylvatica, which suggests that the forest's composition appears to be changing from beech to oak.On the other hand, Benavides et al. (2015) showed that F. sylvatica had an unexpectedly high density of juveniles at lower elevations, suggesting that there may be other reasons for this phenomenon that were not linked to climate change, such as changes of land-use practices and disturbance regimes.With the limited information from our study, we are unable to determine whether the decreased F. crenata regeneration at lower elevations in the Shirakami Mountains may be caused by current global warming or not.To identify the major causes of this pattern, it will be important to examine the effects of other potential factors, such as anthropogenic disturbance.We should also consider aspects of the species' life history, such as changes in flower and seed production (e.g.masting) and in its interactions with other organisms, such as the frequency and severity of damage caused by insects and  seed predators as a function of elevation.Although it is difficult to determine how much the current warming has contributed to this trend, our data agree with previous predictions that the distribution of F. crenata near the low-elevation boundary will shift toward higher elevations as the climate warms (Matsui et al. 2007) and that potential juvenile habitat will be lost at the warmer and snowless sites (Koide et al. 2016).Our results also suggest that if the climate continues to warm in the future, F. crenata regeneration will decrease at lower elevations, where it may be replaced by other species.
In the GLMM analyses, we found a significant positive effect of the F. crenata BA on the number of seedlings (Table 1).This may be because F. crenata trees with a large BA have wide canopies and can therefore produce large numbers of seeds.In the models for saplings, however, the BA of F. crenata had a significant negative effect on the number of saplings in both models, which was the opposite of the result for seedlings.This may be caused by shading by the large mature trees, which would decrease the growth of saplings because F. crenata needs the high light levels that occur in canopy gaps for successful regeneration (Nakashizuka and Figure 3. Relationship between the juvenile-to-canopy tree (J/C) ratio for (a) Fagus crenata, (b) Magnolia obovata, and (c) Acer pictum subsp.mayrii.Note that the y-axis scale differs among the species.Site names are defined in Figure 1 and table S1.S1.
Table 1.Characteristics of the best-fit models for predicting the number of seedlings and saplings of Fagus crenata and the number of saplings of Magnolia obovata and Acer pictum subsp.mayrii in the Shirakami Mountains, Japan, as a function of three explanatory variables: the site's elevation, coverage by dwarf bamboo, and basal area (BA) of individuals with DBH ≥30 cm.The analyses were performed with generalized linear mixed models.The table shows the models with the lowest values of Akaike's information criterion (AIC) and models with values that differ by <2 units from that value; ΔAIC represents the difference in AIC compared to the model with the lowest AIC.The models were constructed based on the results of the tree censuses in 27 plots (each 20 m × 20 m) at the nine study sites.Significance: ***, p < 0.001; **, p < 0.01; *, p < 0.05 (based on t-statistics using Satterthwaite's method).

Response variable
Intercept Numata 1982;Nakashizuka 1987).We hypothesize that the negative effects caused by shading would also be important for the growth of seedlings, but that their numbers would be high because of the abundant seed production by the large trees right after seedling establishment.In the selected models for the number of saplings, the coefficient for the coverage by dwarf bamboo was negative (Table 1).This indicates that dwarf bamboo inhibits the growth of seedlings, which agrees with previous studies (Peters et al. 1992;Abe et al. 2002).
The J/C ratio of M. obovata showed a similar pattern to F. crenata in its relationships between elevation and regeneration, although the relationship was only marginally significant for M. obovata (Figure 3b).This may be because M. obovata has a similar vertical distribution (i.e.200 to 920 m; Matsui et al. 2007) to F. crenata in the Shirakami Mountains, and the young populations in low-elevation areas may also suffer from high mortality caused by current warming or other anthropogenic factors.Our result agrees with Koide et al. (2022), who demonstrated that M. obovata shifted to colder sites by comparing juvenile and adult temperature ranges.However, in the models selected by our GLMM analyses, the M. obovata BA had a significant positive effect on the number of saplings of M. obovata; this effect was the opposite of the relationship for F. crenata.This may be partially due to the sprouting ability of M. obovata, because we often observed young stems of M. obovata originating from the base of large trees.In contrast, the effect of dwarf bamboo was small and not significant in the selected models (Table 1).This may also be due to the high sprouting ability of M. obovata because sprouting stems have a high growth rate and are less sensitive to competition from dwarf bamboo.For A. pictum subsp.mayrii, the coverage of the site by dwarf bamboo and the BA of canopy-size individuals significantly affected its number of saplings, as was the case for F. crenata (Table 1).However, there was no significant relationship with elevation.This may be because its lowest elevation boundary (i.e.0 m according to Matsui et al. (2007)) is much lower than that of our study area.
In the low-elevation areas, the major species that were codominant with F. crenata in the understory layer were A. pictum subsp.mayrii, M. obovata, A. turbinata, and Pterocarya rhoifolia Sieb.et Zucc.(Figure 2).Matsui et al. (2018) predicted that F. crenata is likely to be replaced by Quercus spp., but our data demonstrated that other tree species may also actively regenerate and compete with F. crenata.The species composition under a gap may be different from the ones we found, and additional research will be needed to clarify the mechanisms of species turnover.
The regeneration patterns we obtained in this study were based on 27 plots with closed canopies.We must therefore caution readers that F. crenata regeneration may be more vigorous than would be predicted based on our results if canopy gaps are created by large disturbances (Nakashizuka and Numata 1982;Nakashizuka 1987;Shimano 2006).The J/ C ratios in our study (0.2 to 15.7) were smaller than those (1.2 to 19.4) reported by Shimano (2006) because Shimano focused on well-regenerated areas in canopy gaps.In addition, the regeneration could differ from our results even at the same site, since F. crenata forests have a mosaic structure (Yamamoto et al. 1995;Ida 2000;Henbo et al. 2004).Conversely, it is possible that other species may become even more vigorous than F. crenata in canopy gaps.Considering the mechanisms involved in F. crenata gap dynamics, our plot size may be too small to detect all of this variation following disturbances.

Conclusions
Our results suggest that F. crenata regeneration is poorer at the lower boundary of its vertical distribution than in the middle of its distribution range.If the climate continues to warm in the future, F. crenata regeneration at lower elevations may worsen.It is beyond the scope of our study to examine how much current warming has contributed to the observed regeneration patterns.The decreased regeneration at lower elevations may also be shaped by other factors, such as variation in anthropogenic disturbance or stand history.To more accurately predict the effects of climate change, long-term monitoring studies will be helpful because they provide a better understanding of chronological patterns of turnover and growth rates of the target species.Several monitoring projects are currently underway in the Shirakami Mountains.They are being conducted at a limited number of sites, mainly within the World Natural Heritage site area.We propose that such activities be extended to a wider elevational range, including sites outside the World Natural Heritage site area.Monitoring at low elevations will be especially important because these sites are likely to have been most strongly affected by recent global warming due to the decreased amount of beech regeneration.In addition, these areas are located close to local communities and tourist spots.Participation of citizens in forest monitoring activities would increase public interest in the conservation of beech forests and in the challenges being created by climate change.Such educational activities for F. crenata, a symbol of this World Natural Heritage site, will strengthen the relationships between local people and the ecosystems of the Shirakami Mountains.

Figure 1 .
Figure 1.Map of northern Tohoku, Japan (small map) and map of northern Aomori Prefecture, showing the nine Fagus crenata study sites (large map).The dark and light gray areas in the small map indicate the World Natural Heritage area and the Shirakami Mountains, respectively.The light lines in the large map are municipal boundaries.The base map was created from geographic data provided by the Geospatial information Authority of Japan (https://maps.gsi.go.jp/#8/39.846504/140.789795/&base=blank&ls=blank&disp=1&vs=c0g1j0h0k0l0u0t0z0r0s0m0f1).

Figure 2 .
Figure 2. Stem-diameter distribution of major tree species in forests dominated by Fagus crenata at nine study sites in Shirakami Mountains, Japan.Abbreviations below the x-axis represent the site and elevation range.Figure 1 and table S1 define the abbreviations for the site names.
Figure 2. Stem-diameter distribution of major tree species in forests dominated by Fagus crenata at nine study sites in Shirakami Mountains, Japan.Abbreviations below the x-axis represent the site and elevation range.Figure 1 and table S1 define the abbreviations for the site names.

Figure 4 .
Figure 4. Relationships between elevation and the values of the a (a) and b (b) parameters in the power function constructed based on the DBH class distribution of Fagus crenata: density = a DBH b .Site names are defined in Figure 1 and tableS1.
Elevation Coverage by dwarf bamboo BA of conspecific trees with DBH ≥30 cm AIC ΔAIC