Sustainability of High-Value Timber Species in Mixed Conifer–Broadleaf Forest Managed under Selection System in Northern Japan

Understanding the sustainability of high-value timber species in managed forests provides useful information for the management of these species in the long-run. Using nearly 50 years of census data in long-term permanent plots, we investigated the sustainability of three high-value timber species—monarch birch (Betula maximowicziana Regel), castor aralia (Kalopanax septemlobus (Thunb.) Koidz), and Japanese oak (Quercus crispula Blume)—in cool-temperate mixed forest under a selection system in northern Japan. We used stocking, demographic parameters, and species proportions of these species as measures of sustainability. Results showed that the tree density and basal area of the three high-value timber species increased during the study period. Moreover, the kbasal area increment of these species showed an increasing trend across census periods. However, while no significant differences in the tree mortality of these species were observed, the numbers of in-growth fluctuated across census periods. Increasing trends in species proportions of monarch birch and Japanese oak were observed. Even though there were some fluctuations across census periods, especially in smaller diameter classes, diameter distribution curves of high-value timber species followed a reversed J-shaped pattern. The results revealed that the sustainability measures of high-value timber species can be achieved in forest stands managed under single-tree selection system. In addition, the results also indicated the changing structure and composition of the forest stand. The stocking and basal area increment of conifers decreased while those of broadleaves increased. The proportion of conifers decreased to 33.01% in 2008–2016 from 48.35% in 1968–1978. The results of this study would be useful for adapting silvicultural practices and harvesting practices as well as for simulating various silvicultural and management options for high-value timber species.


Introduction
In uneven-aged mixed conifer-broadleaf forests in northern Japan, monarch birch (Betula maximowicziana Regel), castor aralia (Kalopanax septemlobus (Thunb.) Koidz), and Japanese oak (Quercus crispula Blume) are important producers of high commercial value timber, which is used in the veneer and furniture industries. The supply of high-quality timber from these tree species is exclusively dependent on the cutting of large trees within the mixed forests. Large-sized trees in mixed forests contribute to the structural heterogeneity, dynamics, and functions of the forest ecosystem [1,2], a large fraction of aboveground biomass and carbon storage [3], and play important roles in the rate specific forest management or a single-tree management system, which was recommended in previous studies (e.g., [7,33]).
The aim of this study is, therefore, to assess the sustainability of high-value timber species in mixed forests managed under selection systems. Using 48 years of measurement data, we derived the stocking, demographic parameters, and species proportion of high-value timber species as measures of sustainability. To reach the objective, firstly, we assessed the changes in stocking and demographic characteristics of high-value timber species over time. Secondly, the changes in the species proportion of high-value timber species were assessed. In addition, we also showed how the forest stand structure was changing by assessing the sustainability measures of all conifer and broadleaf species.

Study Site
This study was conducted in the University of Tokyo Hokkaido Forest (UTHF) (43°10-20' N, 142°18-40' E, 190-1459 m a.s.l.), which has an area of 22,717 ha. The UTHF is located in Furano City, central Hokkaido Island in northern Japan (Figure 1). The mean annual temperature and annual precipitation of UTHF are 6.3 °C and 1210 mm, respectively [33]. The ground is covered with snow from late November to early April with a depth of approximately 1 m [34]. Typical vegetation in the UTHF is uneven-aged mixed forests where conifers and broadleaf tree species are the main vegetation cover. The predominant tree species include Sakhalin fir (Abies sachalinensis), Yezo spruce (Picea jezoensis), Japanese linden (Tilia japonica), and painted maple (Acer pictum var. mono) [35]. The forest floor is often occupied by evergreen dwarf bamboo (Sasa senanensis and Sasa kurilensis).

Brief Overview of Single-Tree Selection System at the UTHF
Long-term and large-scale experiments of stand-based forest management system have been conducted since 1958 in the UTHF [35]. Single-tree selection harvest, in which trees are periodically selected and harvested from a large area with cutting cycle of 15 to 20 years and a tree removal rate of 10%-17%, has been implemented as the main silvicultural system in the UTHF [33]. General guidelines for tree marking for selection harvest include tree vitality potential, low quality trees that restrict the growth of seedlings and juveniles, trees with low stem quality, large-sized trees with volume increment that cannot be expected to, or have wood quality that will be deteriorated, or thinning to control stand density, etc. Tree marking operations are conducted by technical staff with long-term operational experiences in tree marking. Attention was paid to be spatially unbiased when

Brief Overview of Single-Tree Selection System at the UTHF
Long-term and large-scale experiments of stand-based forest management system have been conducted since 1958 in the UTHF [35]. Single-tree selection harvest, in which trees are periodically selected and harvested from a large area with cutting cycle of 15 to 20 years and a tree removal rate of 10%-17%, has been implemented as the main silvicultural system in the UTHF [33]. General guidelines for tree marking for selection harvest include tree vitality potential, low quality trees that restrict the growth of seedlings and juveniles, trees with low stem quality, large-sized trees with volume increment that cannot be expected to, or have wood quality that will be deteriorated, or thinning to control stand density, etc. Tree marking operations are conducted by technical staff with long-term operational experiences in tree marking. Attention was paid to be spatially unbiased when marking trees for selection harvest [36]. Larger trees are more likely to be marked than smaller ones, while dominant specie, i.e., A. sachalinensis, is more likely to be selected in tree marking operations [36] even though there is no long-term goal to change the species composition of forest stand. Logging operations typically use felling with chainsaws and skidding of the felled stems by a winch-equipped crawler tractor [27].

Data
Permanent plots were established throughout the UTHF to record long-term growth and stand development for the management of uneven-aged mixed forests [37]. A stratified purposive sampling scheme was employed for the plot establishment to represent the major stand types, soil, and terrain conditions [38]. Harvesting and management of forest stands were monitored using permanent plot data by assessing temporal dynamics of stand structure and fluctuation of forest resources. Within these plots, diameter at breast height (DBH) measurements of all trees with DBH ≥ 5 cm are performed by UTHF staff at regular, in most cases, 5-years interval with 0.1 cm precision. All trees were tagged with identification numbers (IDs) on metal plates nailed to steel rods to ensure that DBH measurements were repeated on the same trees. Plot data included species, DBH, survival status, and harvest of all trees with DBH ≥ 5.0 cm.
In the current study, we used tree census data from 31 permanent sample plots (plot sizes range from 0.22 to 1.00 ha) in the UTHF where the selection system has been implemented. We used the tree census data of the plots measured between 1968 and 2016. Measurement intervals ranged from 3 to 12 years, and the number of measurements in the plots ranged from 8 to 11. The plots were subjected to selection harvest 3 to 5 times during the observation period.

Measures of Sustainability
In order to assess the sustainability of high-value timber species, we considered the stocking and demographic characteristics, and species proportion, as suggested by O'Hara et al. [8]. Stocking included the tree density (number of trees per hectare-N, tree ha −1 ) and basal area per hectare (BA, m 2 ha −1 ). Demographic characteristics included basal area increment (BAI, m 2 ha −1 yr −1 ); number of recruitment or in-growth (N-rec, tree ha −1 ), and number of tree mortality (N-mor, tree ha −1 ). BAI was calculated based on two consecutive measurements of DBH. We excluded negative BAI values in the analysis, as these were assumed to be measurement errors. We also excluded dead trees in the second measurement of two consecutive measurements in BAI calculations. N-rec was calculated as the number of trees that were entered into the minimum 5.0 cm DBH in the second measurement of two consecutive measurements. We considered the dead trees in the second measurement of two consecutive measurements as N-mor, and we counted the number of dead trees between two measurements. In addition, we compared the changes in N, BA, and BAI to understand the general temporal trend of the study forest stand. Changes in the species proportions of high-value timber species as well as the proportion of conifer and broadleaf trees were also analyzed. We also assessed the selection harvest of conifer and broadleaf trees across the census periods.

Data Analysis
For simplicity in the comparison of sustainability measures, we divided the measurement records into 5 census periods; I (1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978), II (1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988), III (1988III ( -1998, IV (1998IV ( -2008, and V (2008-2016). The year of measurement was used to determine the census periods, i.e., all censuses that were measured in the years between 1968 and 1978 were assumed as period I. We tested stand level changes in stocking, demographic characteristics, and proportion of high-value timber species using the Friedman test. For significant Friedman test results (p < 0.05), a post-hoc test was applied to detect which census periods showed significant increase or decrease in stocking, demographic characteristics, or species proportion. A linear mixed effect model was used to quantify whether the stocking, demographic characteristics, and species proportion increased or decreased over time. The model form can be described as the following equation: where Y ij represents the sustainability measures mentioned in Section 2.3 for census year i in plot j, and X ij is census year i in plot j. a 0 and a 1 are fixed effect parameters, and plot j is the random effect parameter for plot j. The symbol ε ij stands for residuals. The significant slope coefficient (p < 0.05) was used as an indicator of either positive or negative trends of the forest characteristics mentioned in Section 2.3 over time [8]. All statistical analysis were carried out in R software version 3.6.1 [39].

Changes in Stocking of High-Value Timber Species
Changes in N and BA are shown in Table 1. Diameter distribution of high-value timber species are shown in Figure 2. During the 48-year period, the N of monarch birch significantly increased (Friedman test, p < 0.05) in the study permanent plots with the mean N of the last census period (2008-2016) being significantly higher than in the first two census periods ( Table 1). The BA of monarch birch also increased, reaching the highest in the last period. Monarch birch BA in the last two periods was significantly higher than in the previous three periods. The N of castor aralia did not increase or decrease significantly (p = 0.75) in any period. However, the BA of castor aralia increased significantly (p < 0.001). The N of Japanese oak did not significantly increase in the first three periods, but it significantly increased in the last two periods. Similar to other two species, the highest BA of Japanese oak occurred in the last period. According to the mixed effect model, the slope coefficients for N and BA for all species showed significant positive trends over time ( Table 2).
where Yij represents the sustainability measures mentioned in Section 2.3 for census year i in plot j, and Xij is census year i in plot j. ao and a1 are fixed effect parameters, and plotj is the random effect parameter for plot j. The symbol εij stands for residuals. The significant slope coefficient (p < 0.05) was used as an indicator of either positive or negative trends of the forest characteristics mentioned in Section 2.3 over time [8]. All statistical analysis were carried out in R software version 3.6.1 [39].

Changes in Stocking of High-Value Timber Species
Changes in N and BA are shown in Table 1. Diameter distribution of high-value timber species are shown in Figure 2. During the 48-year period, the N of monarch birch significantly increased (Friedman test, p < 0.05) in the study permanent plots with the mean N of the last census period (2008-2016) being significantly higher than in the first two census periods ( Table 1). The BA of monarch birch also increased, reaching the highest in the last period. Monarch birch BA in the last two periods was significantly higher than in the previous three periods. The N of castor aralia did not increase or decrease significantly (p = 0.75) in any period. However, the BA of castor aralia increased significantly (p < 0.001). The N of Japanese oak did not significantly increase in the first three periods, but it significantly increased in the last two periods. Similar to other two species, the highest BA of Japanese oak occurred in the last period. According to the mixed effect model, the slope coefficients for N and BA for all species showed significant positive trends over time ( Table 2).    For the purpose of comparison, Figure 3 shows changes in the N and BA of forest stand, conifer, and broadleaf across census periods. The Friedman test showed significant differences in both N and BA over time (p < 0.05). In addition, Figure 3 also shows the changes in the N and BA of conifer and broadleaf species during the 48-year period. The figure shows increasing trends in broadleaf N and BA, while decreasing trends can be observed for conifer N and BA. Friedman test also showed significant differences (p < 0.05). The results of the mixed effect model also confirmed the significant and positive trends in N and BA over time. However, significant and negative trends were observed for conifer N and BA, while a positive trend was observed for broadleaf species (Table 3).  For the purpose of comparison, Figure 3 shows changes in the N and BA of forest stand, conifer, and broadleaf across census periods. The Friedman test showed significant differences in both N and BA over time (p < 0.05). In addition, Figure 3 also shows the changes in the N and BA of conifer and broadleaf species during the 48-year period. The figure shows increasing trends in broadleaf N and BA, while decreasing trends can be observed for conifer N and BA. Friedman test also showed significant differences (p < 0.05). The results of the mixed effect model also confirmed the significant and positive trends in N and BA over time. However, significant and negative trends were observed for conifer N and BA, while a positive trend was observed for broadleaf species (Table 3).

Changes in the Demographic Characteristics of High-Value Timber Species
Changes in the demographic parameters of high-value timber species are shown in Table 4. The BAI of monarch birch in the last two census periods was significantly higher than the first three periods (p < 0.001). A significant higher BAI of Japanese oak (p < 0.001) was found in the last two measurement periods with the highest increment in the fourth census period. An increasing trend in the BAI of castor aralia was also observed. The mixed effect model also revealed a positive slope

Changes in the Demographic Characteristics of High-Value Timber Species
Changes in the demographic parameters of high-value timber species are shown in Table 4. The BAI of monarch birch in the last two census periods was significantly higher than the first three periods (p < 0.001). A significant higher BAI of Japanese oak (p < 0.001) was found in the last two Forests 2020, 11, 484 7 of 13 measurement periods with the highest increment in the fourth census period. An increasing trend in the BAI of castor aralia was also observed. The mixed effect model also revealed a positive slope coefficient for the BAI of all species ( Table 2). The Friedman test showed no statistically significant difference between census periods for N-mor (Table 4). However, a significant positive slope was observed for monarch birch and Japanese oak, meaning that the N-mor of monarch birch and Japanese oak may be likely to increase over time ( Table 2). In all target species, an increased N-rec was observed in the second, third, and fourth census periods. However, a significant decline in the N-rec was observed in the last census period (Table 4), and no significant positive or negative slope coefficients were observed for any species in the results of the mixed effect model ( Table 2).  Figure 4 shows changes in the BAI, N-mor and N-rec of all species. The total BAI of all species showed no statistically significant difference in the first three census periods. However, it was significantly higher in the last two periods than in the first three periods (p < 0.05). The total BAI of broadleaf species increased while the conifer BAI decreased in the last census period. These results are also confirmed by the results of the mixed effect model in which a significant negative slope was observed for conifer and a positive slope was observed for broadleaf. The broadleaf N-rec that entered into the 5.0 cm DBH was higher than that of the conifer N-rec. For both conifer and broadleaf species, the N-rec decreased from the third census period (Figure 4). The number of total N-mor significantly increased since third census period. When comparing the conifer and broadleaf N-mor, the broadleaf N-mor increased from the third census period. However, there was no significant increase or decrease in conifer N-mor after the third census period. The mixed effect model also revealed an increase in the total N-mor over time. In addition, the broadleaf N-mor also increased over time. Table 5 shows the changes in species proportion of high-value timber species. Significant differences in the proportion of Japanese oak were observed (p < 0.001), while no significant differences were observed for monarch birch (p < 0.22) or castor aralia (p < 0.31) across census periods. However, a significant positive trend in the monarch birch proportion was observed over time (Table 2). Similarly, Table 2 shows a significant positive trend in the Japanese oak proportion. Table 5 also shows the changes in the proportion of conifer and broadleaf across census periods. The proportion of conifer declined after the first census period while that of broadleaf increased after the first period. These trends were also confirmed by the positive coefficients of mixed-effect models (Table 2). species, the N-rec decreased from the third census period (Figure 4). The number of total N-mor significantly increased since third census period. When comparing the conifer and broadleaf N-mor, the broadleaf N-mor increased from the third census period. However, there was no significant increase or decrease in conifer N-mor after the third census period. The mixed effect model also revealed an increase in the total N-mor over time. In addition, the broadleaf N-mor also increased over time.     Table 6 shows the total N and BA harvest of high-value timber species. The small N and BA of high-value timber species were harvested during the 48-year period. Table 7 shows the total N and BA harvest of conifer and broadleaf across census periods. The number of broadleaf trees harvested in the first census period was higher than the number of conifer trees harvested. Starting from the second census period, more conifer trees than broadleaf trees were harvested. According to Table 6, the BA harvest for conifer was larger than the BA harvest for broadleaf in all periods. The BA harvest for broadleaf in the first census period was significantly higher than in the later four periods, with the lowest BA harvest occurring in the third period. The largest BA harvest for conifer occurred in the fourth census period.

Discussion
In this study, we investigated the long-term changes in stocking, demographic characteristics, and species proportion of high-value timber species as measures of sustainability. We used these parameters as criteria and indicators to evaluate the sustainability of high-value timber species as proposed by O'Hara et al. [8]. However, they used these criteria and indicators to evaluate the stand level sustainability of even-aged and uneven-aged forest management. The changes or consistency of these parameters would be useful for establishing sustainable forest management by adjusting the tree marking for harvesting.
The main common characteristics of the targeted high-value timber species in this study were increases in N and BA. The increasing trend was also observed in the total N and total BA of forest stands. An increase in the N of forest stands managed under the selection system has occurred in the last few decades in other parts of the world with different environment and forest types. For example, Klopcic et al. [20] found an increased total number of trees in their study in Slovenia. Compared with an unmanaged stand, a higher mean density and basal area in the managed stand was also reported by Young et al. [23] in USA. Moreover, a decreasing tree density and basal area in an unmanaged stand was reported by Ediriweera et al. [24] in their study in mixed-dipterocarp forests over a 40-year period.
One possible reason for increasing the N and BA of forest, including those of high-value timber species, would be due to some major disturbance, including natural (e.g., strong typhoon) and anthropogenic (e.g., selection harvest) factors [26]. In the mixed conifer-broadleaf forest in northern Japan, a large typhoon occurred, causing widespread canopy opening [40] in some plots. Selection harvesting was carried out three to four times during the 48 years period. The increasing trend of both total N and BA of forest stand was greatly contributed by broadleaf species. These increasing trends in both N and BA are consistent with previous studies in mixed conifer-broadleaf forests in northern Japan such as that done by Yoshida et al. [18]. Those studies also highlighted the increasing trend of broadleaf tree density in mixed conifer-broadleaf forests managed under selection system.
Reversed J-shaped diameter distributions were observed for all target high-value timber species (Figure 2). Diameter distribution curves indicated an increasing number of smaller diameter class trees for all species across census periods. Owari et al. [35] also reported a large number of smaller diameter class trees in northern Japanese mixed conifer-broadleaf forests. More fluctuation in size structure of monarch birch and Japanese oak were observed.
Similar to N and BA, the mean BAI of high-value timber increased over time. The increasing trend was also observed for the mean stand BAI even though it was not significantly different in the first three census periods. The results of regression modeling suggested temporal trends of BAI for high-value timber species and the total BAI of forest stands. In terms of the stand level BAI, a significant positive trend was observed. In northern mixed conifer-broadleaf forests, similar results have been reported for the broadleaf species (e.g., [18,41]).
The mean N-mor across census periods showed no significant difference for all high-value timber species. However, Japanese oak N-mor increased over time. Furthermore, no significant increasing or decreasing trends in N-rec were observed for high-value timber species across census periods in this study. These results are also consistent with previous studies by Hiura et al. [41]. For all target species, a significantly lower N-rec was observed in the first and last census periods, while a higher N-rec was observed in the second, third and fourth census periods. This trend can also be observed for the total N-rec of all species. This pattern may also be explained by a large disturbance caused by a large typhoon in 1981 (second census period) which created a canopy opening causing high light availability [42] in some plots, favoring natural regeneration. As a result, higher N-rec in these plots would be expected following a large natural disturbance. This was also highlighted by previous studies in the mixed conifer-broadleaf forest in Northern Japan. Yoshida et al. [18], for example, reported that a higher N-rec of castor aralia would be expected after a disturbance, and harvesting treatment in their study. In addition, Takahashi et al. [43] reported that Japanese oak may tend to regenerate after a large disturbance before the establishment of other species.
The results in this study also clearly show another important temporal trend of the forest stand. The forest composition and structure in the study area changed over time after the first census period (1968 to 1978). The proportion of conifer in the first census period was 48.35%, and it decreased to 33% in the last census period. An increasing broadleaf proportion might contribute to an increased N of high-value timber species (Table 2). Broadleaf N and BA exceeded those of conifer after the first census period, wherein the N and BA of conifer were larger than those of broadleaf. Even though there might be several reasons for this, one of the possible reasons for decreased conifer N and BA would be due to a bias in tree marking for harvesting. Table 6 showed that the total harvested basal area for conifer always exceeds that of broadleaf in all census periods. Owari et al. [36] indicated that even though spatially unbiased tree marking was obtained, tree species that were marked for selection harvest in our study area was mostly conifer species. In addition, they observed that the trees marked for harvesting were larger than the unmarked trees.
In terms of BAI, conifer showed a decreasing trend over time (Table 3). Even though the mean BAI of conifer exceeded the mean BAI of broadleaf throughout the census periods, both were almost similar in the last census period (0.379 m 2 /ha for conifer and 0.369 m 2 /ha for broadleaf). In addition, the mean BAI of conifer was significantly lower in the last census period than in the first four periods. Yoshida et al. [18] reported that the growth of shade-intolerant broadleaf species was less than that of most common conifer species (i.e., A. sachalinensis), even after selection harvesting. However, our study found that even though a larger mean BAI was observed for all conifer species, its BAI showed a decreasing trend while the total BAI of broadleaf species exhibited increasing trend. Harvesting of more conifer trees may contribute lower mean value in total BAI. Our results are in line with a previous study by Hiura et al. [41] in mixed conifer-broadleaf forest in northern Japan. According to this trend, the future forest composition of the study area would be more of a broadleaf-dominating type in terms of both N and BA. A similar decreasing trend in the conifer proportion has been reported in different regions [20,21].
The results of the sustainability measures in this study revealed that there have been inconstancies in these measures over time. However, such inconsistencies could relate to a number of reasons. A recent study [41] in mixed conifer-broadleaf forest with very little human disturbance in northern Japan revealed that changing climate conditions such as an increased temperature, precipitation, and decreased snowfall and snow cover period have led a to reduction in growth rate of conifer and an increasing in that of broadleaf species. Moreover, these inconsistencies in sustainability measures due to a changing climate have been widely reported in different regions [44,45]. Similar to the results by O'Hara et al. [8], the results of this study indicate that a single-tree selection system is more of a dynamic entity.
The sustainability measures described in this study would be useful for adjusting forest management activities, and various silvicultural activities, which could lead to consistency in sustainability measures in different forest types. Through the understanding of sustainability measures used in this study, forest management can maintain the stocking of uneven-aged forest stand over time, BAI can be balanced by tree removals, and recruitment can be assessed whether it is sufficient. In addition, it would provide information for forest management operations such as stocking control, which is central to uneven-aged silviculture. Many stocking control approaches have been developed including reversed J-shape diameter distribution, selection system or plenter system, stand density index, and leaf area allocation, etc. [46]. Sustainability measures can be achieved through these stocking control approaches by removing those trees that surpass or maintain those of limited numbers in a certain diameter class.

Conclusions
In this study, we examined the sustainability of high-value timber species in mixed conifer-broadleaf forest managed under selection system. Changes in the stocking, demographic characteristics, and species proportion of high-value timber species over a 48-year period have been used as measures of sustainability. These measures could provide useful information for their management in the long-run. In addition, we examined the general temporal trend of a selection forest stand. The main common characteristics of high-value timber species were increases in the tree density, basal area, and BAI across the census period. In terms of tree mortality, no significant differences were observed among census periods with no significant downward or upward trends. High fluctuation in the number of in-growth also occurred. Other important long-term characteristics of the forest stands are the changes in forest structure and composition to broadleaf forest in terms of the tree density and basal area. Even though some fluctuations in sustainability measures were observed, the results indicated that sustainability of high-value timber species was achieved under a single-tree selection system. The results of this study would be useful for adapting silvicultural practices and harvesting practices such as single-tree selection, and also to maintain desired stocking of the forest stand. For example, attention should be paid to possible changes in species proportion and diameter distribution of remaining forest stand when marking trees for harvesting. The influence of natural and anthropogenic factors on long-term changes in stocking and demography of high-value timber species should be analyzed further.