Intra-specific variation in mortality of even-aged Cryptomeria japonica (L. f.) D. Don. forests can be explained using relationships among long-term stand characteristics

Context

Understanding tree mortality is critical for sustainable forest management. Long-term tree mortality may differ depending on the stand development process and can be influenced by forest management regimes. Logistic regression is widely used to explain tree mortality based on site productivity, age, size, and competition. However, the explanatory variables are interrelated. Thus, we attempted to explain long-term Japanese cedar tree mortality by considering interrelated variables.

Aims

The aim of this study was to elucidate the direct and indirect effects of site productivity, age, individual size, and competition on the long-term mortality of Japanese cedars.

Methods

Data were collected from 5130 even-aged Japanese cedar trees over approximately 50 years. We compared each variable between dead and living trees. We then constructed a mortality model using a conventional logistic approach and selected the best model for the stepwise methods. Finally, we applied a piecewise structural equation model (SEM) to identify these variables’ direct and indirect effects. We compared the conventional logistic model and piecewise SEM models and discussed the advantage of applying the SEM models.

Results

Annual mortality was approximately 4% in the most fertile stands, increasing gradually with decreasing site fertility. Dead tree size and competition status differed according to age and site productivity. Competition, individual size, and stand density were selected for the best logistic model (area under the curve (AUC) = 0.74, Brier score = 0.042), whereas age and site productivity were not (p > 0.05). The piecewise SEM results showed that age and site productivity indirectly affected tree mortality through individual size and stand density (Fisher’s C = 4.569, p = 0.102).

Conclusion

Long-term Japanese cedar tree mortality can be explained by individual size and competition as direct influencing factors and age and site productivity as indirect influencing factors. This indicated that hidden factors cannot be explained using the conventional logistic approach. Further studies are required to explore the potential factors contributing to tree mortality thoroughly.

Individual size and competition directly affected long-term Japanese cedar (Cryptomeria japonica (L. f.) D. Don.) mortality, whereas age and site productivity indirectly affected the mortality through individual size and competition. This result suggested that there are potential factors that could not be explained using a conventional logistic approach.

Understanding tree mortality is critical for sustainable forest management. Tree mortality affects forest yield, structure, dynamics, habitat, and biodiversity (Tilman 1994). It is influenced by various factors, such as competition (Weiskittel et al. 2011), disease (Brasier, 1996), pathogens (Breece et al. 2008), fire (Stephens and Moghaddas 2005), and wind damage (Nagel and Diaci 2006). Tree mortality is classified as regular or irregular based on the cause of death (Lee 1971). Regular mortality is caused by competition, whereas irregular mortality is caused by disease, pathogens, and wind damage (Lee 1971). As irregular mortality results from sudden disturbances, it is difficult to predict. Most previous studies have focused on regular mortality when modeling tree mortality (Weiskittel et al. 2011), a focus that was also adopted in this study.

Tree longevity differs among species, as observed, for example, in pioneer, climax, and wood species, which have different life strategies (e.g., shade tolerance and growth speed) and also vary depending on sites and silvicultural treatments, even for the same species. Old trees are often found in harsh environments and poor sites, such as high elevations or extremely dry areas (Loehle 1988; Das et al. 2016). Consequently, old trees may be smaller than other trees of the same species, or, as expressed by age, the largest trees are not the oldest (Piovesan and Biondi 2021). However, the social or crown class, defined by crown width, tree vigor, and vertical stand structure, is widely adapted to even-aged pure forests (Burkhart et al. 2019). In general, suppressed trees in even-aged stands are smaller and have a high mortality risk because their growth is limited by competition with other trees, eventually leading to death (Franklin et al. 1987). This implies that even-aged trees die differently in their surrounding environments.

Numerous studies have used tree or stand mortality models in even-aged stands to describe the underlying biological reality. Binomial logistic regression is widely used to predict mortality and analyze various factors. (Etzold et al. 2019; Lee and Choi 2019; Ningre et al. 2019; Zhang et al. 2020; Pretzsch et al. 2022). This suggested that tree mortality rates tend to be higher when site productivity is high, trees are young and small, and competition is high (Chen et al. 2008; Kuehne et al. 2022). The logistic regression model directly incorporated these variables as independent mortality factors (Cao 2017; Saud et al. 2022). However, age affects individual size; that is, tree size increases with age (Thomas 2013), and site productivity affects tree size; that is, trees grow faster on sites with high site productivity (Nishizono 2010). This indicates that the factors that potentially affect tree mortality are interrelated. The logistic approach could predict tree mortality; however, this approach, which considers only the direct effect, might be insufficient for understanding real tree mortality. A different approach that considers the relationships between factors is needed to organize the evidence on the processes underlying mortality. Piecewise structural equation model (SEM) distinguishes direct and indirect factors and applies a logistic model (Lefcheck 2016). It is expected that this method will not only be able to build a prediction model for mortality but also to clarify the relationship between the explanatory variables, leading to a deeper understanding.

The Japanese cedar (Cryptomeria japonica (L.f.) D. Don.) is one of the most important species used in Japan’s forestry industry. It is characterized by straight stems, and the stands are typically planted at the same time. Many studies on even-aged Japanese cedar have focused on tree growth (Nishizono et al. 2008; Kitagawa et al. 2018). Although these factors are known to affect tree mortality over limited periods, their effects on long-term tree mortality have not yet been clarified (e.g., Osumi et al. 2000; Nishizono et al. 2008; Fukumoto et al. 2022). An adequate description of long-term tree mortality can assist efforts toward efficient wood productivity and sustainable forest management.

Our main objective was determining how interrelated variables, such as site productivity, age, individual size, and competition, influence long-term Japanese cedar mortality. First, we used the collected data to compare the variables of dead and living trees. We then constructed a mortality model using a logistic approach to identify the variables that directly affected tree mortality. Finally, piecewise SEM was used to identify their relationships with direct or indirect effects on tree mortality. From these results, we clarified (1) how site productivity, age, individual size, and competition affect tree mortality and (2) the advantage of applying the logistic model through the SEM model rather than a simple logistic model.

Locations of permanent sampling sites used for investigating long-term cedar mortality in western Japan. In total, 16 permanent plots were established. Site A: Asagihara, Site B: Nishimatahigashimata-yama, Site C: Nakanokawa-yama, and Site D: Kudarukawa-yama. The number indicates plot No (c.f. Table 1). The solid line indicates a difference in initial planting density, and the dashed line indicates the plots established at different slopes

Full size image

2.1 Data collection

The data were collected from four study sites located in the Shikoku region, western Japan (Fig. 1), where the mean annual temperature is 16.2 °C and the mean precipitation is 1322.5 mm/year. The study sites were national forests located in Asagihara (Site A), Nishimatahigashimata-yama (Site B), Nakanokawa-yama (Site C), and Kudarukawa-yama (Site D) (Fig. 1, Table 1 ) (Fukumoto et al. 2022). The soil type in all sampling plots was brown forest soil. These sites were originally established to study the effects of thinning and initial planting density on the growth of Japanese cedar. The trees were planted at each site between 1950 and 1964. Bare root seedlings were planted individually on hill slopes using a hoe, and the same tillage treatment was conducted in all plots. Two to six study plots with areas ranging from 350 to 2270 m² were established between 1960 and 1972. In total, there were 16 study plots in the study sites (Fig. 1 and Table 1 ). Stand densities ranging from 1308 trees/ha to 5250 trees/ha were established at each study site. There were seven high-density plots (< 3000 trees/ha) and nine low-density plots (> 3000 trees/ha) during the first measurement. The first measurements were obtained when the plots were 10–27 years old. A census was carried out approximately every five years, and the measurements were repeated five to ten times. The site index (SI) values for each plot ranged from 15.6 to 27.7 m. In this study, the SI was calculated by measuring the upper mean tree height (250 trees/ha) of 40-year-old trees in each plot (West 2015). There were four high SI plots (SI > 25), six middle SI plots (20 < SI < 25), and six low SI plots (SI < 20).

The study sites had thinned and unthinned plots, with the former being thinned to suppress trees a maximum of two times over approximately 60 years. Residual tree diameter at breast height (DBH) and height were measured in each plot during each census. DBH was measured for all trees in each plot, with the breast height located 1.2 m above the ground and marked with a line. The heights of approximately 30 randomly selected trees were measured in each plot, and the heights of the remaining trees were estimated using the Näslund equation (Nigul et al. 2021). The Blum-Leisse instrument was used to measure height before 2000, whereas the Vertex III (Haglöf Sweden AB, Långsele, Sweden) was used thereafter. The dead and thinned trees were recorded during each field measurement period. Tree mortality was represented by either 0 (dead trees) or 1 (living trees). The measurement data of dead trees used previously recorded data because the size of the dead trees was not measured. In total, we collected data from 5130 trees.

2.2 Measurement of each influencing factor for both dead and living trees

We obtained the mean values of each factor for dead and living trees at sites with different productivity levels and age classes to clarify the effects of site productivity, age, size, and competition on dead and living trees. We then defined site productivity levels as high (SI > 25), middle (20 < SI < 25), and low (SI < 20). Similarly, age was classified as older (age > 50), mature (30 < age < 50), or young (age < 30). First, the mean mortality of the different classes was calculated. We also calculated the mean DBH and competition index for living and dead trees in the different classes and then applied a t-test to compare each factor between dead and living trees.

2.3 Individual mortality model

We applied a generalized linear mixed-effects model with a logistic function to our dataset to model the mortality rate of individual trees. Because the intervals differed by varying degrees between the censuses, we used an Exposure method based on Shaffer (2004). The exposure model is a modified logistic regression incorporating time interval t into the link function. The details of this method are provided in Shaffer (2004). The probability mortality rate \({M}_{I,i,j}\) of the subject tree i in the Ith plot of the jth measurement is expressed as

where \({Age}_{I,j}\) is the stand age of the Ith plot of the jth measurement; \({DBH}_{I,i,j}\) is the DBH value of subject tree i in the Ith plot of the jth measurement; \({SI}_{I}\) is the SI at the Ith plot;\({D}_{I,j}\) is the stand density of the Ith plot in the j ^th measurement; and \({BAL}_{I,i,j}\) (basal area of trees larger than the subject tree) represents the distance-independent competition indexes. BAL is expressed as:

where \(BAL_{I,i,j}\) is the sum of the basal area of the trees competing against subject tree i in the Ith plot of the jth measurement; \({DBH}_{c}^{2}\) is the sum of the basal area of the competition trees with DBH values higher than those of the subject tree i in the Ith plot of the jth measurement (Tenzin et al. 2017); \(a_1-a_5\) are individual parameters; and \({\varphi }_{a,I}\) is a normally distributed random parameter for the Ith plot. \({Age}_{I,j}\), \({SI}_{I}\), and \({D}_{I,j}\) were stand-level factors (n = 16), whereas \({BAL}_{I,i,j}\) and \({DBH}_{I,i,j}\) were tree-level factors (n = 5130). We then used the data for all dead and living trees for statistical analysis and estimated each parameter using the lme4 package in R v.4.1.3. (Bates et al. 2015; R Core Team 2022).

Akaike’s information criterion (AIC) and the area under the curve (AUC) were used to select explanatory variables and evaluate model performance, respectively (Godeau et al. 2020; Hanley and McNeil 1982). The model with the lowest AIC was defined as the best model. The coefficient values were used to evaluate whether the explanatory variables positively or negatively affected tree mortality. AUC values were calculated by drawing a receiver operating characteristic curve (Pencina et al. 2008) and ranged from 0 to 1, with 1.0 indicating perfect distinction. In this study, we calculated AUC values using only fixed-effect variables. Similarly, we calculated the Brier score (Brier 1950), which approached zero as the model accuracy increased. We used the pROC (Robin et al. 2011) and DescTools (Signorell 2020) software packages to calculate the AUC and Brier scores, respectively.

2.4 Path analysis

SEM allows to model multivariate relationships by combining two or more structural models. This approach is widely used to quantify ecological systems and is a useful tool for quantifying the direct and indirect effects on target variables. Path analysis was used to assess the direct and indirect relationships between predictors and tree mortality. We then adapted piecewise SEM using the R package for piecewise SEM (Lefcheck 2016). Piecewise SEM is based on the traditional SEM and can incorporate a wide model structure, distribution, and assumptions. We constructed piecewise SEMs based on the pre-analysis, and the best model was selected using AIC, p-values, and Fisher’s C statistics (Lefcheck 2016). The models were as follows:

The response variables included DBH, stand density, and mortality. Individual size (\({DBH}_{I,i,j}\)) and stand density (\({D}_{I,j}\)) were selected as the key response variables for mortality (\({M}_{I,i,j + 1})\). \({\varphi }_{c,I}\), \({\varphi }_{d,I}\) and \({\varphi }_{e,I}\) are normally distributed random parameters for the I ^th plot. We assumed that (1) DBH is positively affected by age and negatively affected by BAL and stand density; (2) stand density is negatively affected by age, BAL, and DBH; and (3) mortality is positively affected by BAL, SI, and stand density, and negatively affected by age and DBH. Equations (4) and (5) follow a normal distribution, whereas Eq. (6) follows a binomial distribution with Exposure methods. Equations (4) and (5) were constructed using the linear mixed model via the nlme package (Pinheiro et al. 2021), whereas Eq. (6) was constructed using glmm [generalized linear mixed model] in the lme4 package (Bates et al. 2015). All response variables were standardized using mean and variance. Next, we calculated the effect of each predictor on the response variable using p-values and Fisher’s C statistics (Lefcheck 2016), with p > 0.05 and a small Fisher’s C indicating a good model fit. We also calculated the marginal and conditional R ² values to evaluate the accuracy of Eqs. (4) and (5) (Nakagawa and Schielzeth 2013). When piecewise SEM detects variables that indirectly affect tree mortality, it indicates that some variables cannot be represented using the usual logistic approach. In this case, we concluded that the SEM has an advantage over conventional logistic models.

3.1 Effect of each factor on mortality

The mean mortality rate was the highest in young trees in the high SI class (4.5%), and the mean mortality of mature and older plots was 3.3% Fig. 2 (a). The mean mortality rate was highest in mature trees (5.8%) and lowest in young trees (0.2%) in the middle SI class. The mean mortality rate was lowest in young trees (1.3%), whereas the rates for mature and older trees were 3.9% and 2.7%, respectively, in the low SI class. The mean mortality rate differed widely with the SI in young trees, whereas it remained constant in older trees, regardless of the SI.

There was a significant difference in mean BAL between living and dead trees in all age classes in the high SI class (p < 0.001; Fig. 2(b)), and the mean BAL was larger in dead trees (8.1–11.8) than in living trees (4.1–7.2). There was a large difference in the mean BAL between living and dead trees in the high SI class, regardless of age. In the middle and low SI classes, there was no significant difference in this parameter between the living and dead trees at maturity (p > 0.05). In the middle SI class, the mean BAL of living trees ranged from 4.9 to 10.0, and that of dead trees ranged from 6.1 to 13.7, while in the low SI class, it ranged from 1.8 to 7.1 in living trees and from 2.9 to 8.5 in dead trees.

The mean DBH of dead trees was significantly lower than that of living trees in all SI and age classes (p < 0.001; Fig. 2 (c)). It ranged from 10.4 to 22.4 cm in the high SI class and from 9.6 to 16.9 cm in the middle SI class, while it ranged from 3.6 to 9.7 cm in the low SI class and was smaller than that in the other SI classes. In addition, it reached 10 cm in younger ages in the high and middle SI classes and in older ages in the low SI class.

Comparison of the variables in all study plots and field surveys for different site index (SI) and stand age classes; (a) mean mortality rate, (b) basal area of trees larger than the subject tree (BAL), (c) diameter at breast height (DBH). SI classes were defined as high (SI > 25), middle (20 < SI < 25), and low (SI < 20). Age classes were defined as older (age > 50), mature (30 < age < 50), and young (age < 30). The error bar indicates standard error. **** indicates a significant difference (t-test, p < 0.001). The error bar indicates the standard deviation

Full size image

3.2 Best model for estimating long-term mortality in Japanese cedar with the conventional logistic model

BAL, DBH, and D were selected using a stepwise AIC, whereas age and SI were excluded (Table 2). The coefficients of the three parameters were 1.210, 2.471, and 1.549, and their SE ranged from 0.113 to 0.158. BAL had a positive effect on long-term mortality, whereas DBH and D had negative effects. By contrast, the SE values for age and SI were 0.165 and 0.598, respectively. The AUC and the Brier score were 0.74 and 0.042, respectively.

3.3 Effect of predictors on tree mortality with the piecewise SEM

The piecewise SEM results showed that the models were well-fitted (Fisher’s C = 4.569, p = 0.102). First, DBH was affected by age, BAL, and SI, with SE values ranging from 0.0023 to 0.0935 (Table 3 and Fig. 3). The coefficient values of these three parameters were 0.974, 0.751, and 0.846, respectively. Additionally, the marginal R ² and conditional R ² were 0.79 and 0.90, respectively. D was negatively affected by age, BAL, and DBH (marginal R ² = 0.09, conditional R ² = 0.91), with SE ranging from 0.0036 to 0.0040. The coefficient values of these three parameters were − 0.191, − 0.016, and − 0.095, respectively. Finally, the mortality rate was affected by BAL, DBH, and D (cf. the previous paragraph). BAL had a positive effect, whereas DBH and D negatively affected the mortality rate. Age and SI indirectly affected tree mortality rates through DBH and D (Fig. 3).

Direct or indirect relationship between tree mortality and each of the following factors: SI, age, BAL, DBH, and stand density. Arrows indicate unidirectional effects among variables. The red and black arrows indicate positive and negative effects, respectively. The gray-dotted line indicates non-significant effects (p > 0.05). The values in the figure show the coefficients of each parameter based on the piecewise structural equation model (Table 2)

Full size image

At sites with high productivity, the mortality rate was higher in the younger age group than in the other SI classes (Fig. 2 (a)). A large difference in the basal area of trees larger than the subject tree (BAL) was observed between living and dead trees, indicating high competition in the stands (Fig. 2 (b) ). Conversely, the mortality rate gradually increased with age at sites with low productivity. This suggested that stand development is slower in sites with low than with high productivity, which is consistent with the results of a Pretzsch et al. (2023). Moreover, BAL and DBH values showed that the suppressed trees had died (Fig. 2 (b)and (c) ). One factor identified as having the highest mortality risk is limited light (Zuleta et al. 2022), thus small trees are at high risk. These results suggested that tree death occurred at different ages and under different competition conditions, even though the DBH remained the same.

The logistic model results showed that tree mortality rate was affected by competition and tree size but not by age and site productivity (Table 2). These results indicate that Japanese cedar tree mortality can be explained by individual size and competition with conventional logistic approaches. Competition had a positive effect and individual size negatively affected tree mortality. This indicated that tree mortality tended to increase with more intense competition and smaller tree sizes. This result supported the general trend (e.g., Saud et al. 2022). Our results also showed that DBH was the most critical variable for mortality because its coefficient was the largest in Eq. 7 (Table 3). Takata and Kobayashi (1983) indicated that individual tree size plays a more important role in mortality than does intraspecific competition in Japanese cedar growth. Many studies have also indicated that individual size is an expected high interpretability factor for the mortality model; therefore, sometimes tree age has been replaced by individual size when stand or tree age is not available, such as in mixed or uneven-aged stands (e.g., Etzold et al. 2019). Moreover, tree mortality was lower with high stand density, and this result showed a different trend in other studies (e.g., Yang et al. 2024). Loehle (1988) showed that tree growth declined to extend longevity at sites with limited water and nutrient availability. Based on this theory, Chao et al. (2008) indicated that the negative correlation between stand density and mortality reflects a fundamental trade-off between resource acquisition and investment survival. In this study, tree mortality was lower at sites with low productivity and high stand density. Therefore, stand density may have negatively affected tree mortality in the logistic model. In contrast, age and site productivity were not significant variables in our logistic model (Table 2). This indicated that age and site productivity are insufficient for directly describing the long-term mortality of Japanese cedar trees. Many studies have directly inserted age and site productivity as explanatory variables for the logistic mortality model (e.g., Zhang et al. 2020), and these were significant variables for their models (Timilsina and Staudhammer 2012; Cruickshank 2017; Cao 2019). However, individual size and competition status differed by age and site productivity (Fig. 2), and the usual logistic approach might not sufficiently explain tree mortality under different stand development phases.

The piecewise SEM results showed that tree mortality rate was explained by the direct effects of tree size, competition, and stand density and the indirect effects of site productivity and age (Fig. 3). This indicated that tree mortality increases with smaller tree sizes and intense competition; however, tree size and competition are affected by site productivity and age. Pretzsch et al. (2022) indicated that heterogeneity in size structure increases in poor sites because of high competition and small tree death at rich sites. Loehle (1988) pointed out that the relative growth rate of trees decreases to expand their life span in poor sites, and the ratio of photosynthesis to respiration decreases due to increased respiration demands from the supporting tissue. These results suggested that piecewise SEM models adequately reflect the different situations in which tree mortality occurs at different site productivities and ages. Conventional logistics approaches overlook the indirect factors of site productivity and age. If we analyzed tree mortality using data collected from long-term, varied locations and silvicultural treatments, the piecewise SEM model would have an advantage in understanding tree mortality. Piecewise SEM has been widely used for factor analyses recently (Stenegren et al. 2017; Bomfim et al. 2021). This approach can identify the relationships among interrelated variables (Lefcheck 2016) and can potentially be applied to various factor analyses. This study focused on even-aged pure stands. Similar results may not be obtained for other forest types, such as mixed- or uneven-aged forests. Therefore, future studies should apply these models to various forest types.

This study clarified that the interrelated variables of site productivity, age, individual size, and competition influence long-term Japanese cedar mortality. The collected data showed that dead tree size and competition status varied depending on age and site productivity. However, age and site productivity were not significant variables in the conventional logistic regression model. They indirectly influenced long-term Japanese cedar tree mortality through individual size and competition. Our findings highlighted the fact that important variables are overlooked in the conventional logistic approach. Therefore, piecewise SEM has an advantage when adapted to factor analysis. The explanatory variables corresponding to the data should be carefully selected when describing and understanding tree mortality.

The datasets generated during and/or analyzed during the current study are publicly available from the corresponding author at reasonable request.

We thank the editors and two anonymous reviewers for their constructive comments.

Code availability

Not applicable.

This work was supported by research grant #201804 of the Forestry and Forest Products Research Institute.

Authors and Affiliations

1. Kyushu Research Center, Forestry and Forest Products Research Institute, 4-11-16 Kurokami, Chuo, Kumamoto, Kumamoto, 860-0862, Japan

   Keiko Fukumoto

2. Forestry and Forest Products Research Institute, 1 Matsunosato, Tsukuba, Ibaraki, 305-8687, Japan

   Tomohiro Nishizono & Fumiaki Kitahara

Contributions

Conceptualization: Keiko Fukumoto; Methodology; Keiko Fukumoto; formal analysis and investigation: Keiko Fukumoto, Tomohiro Nishizono, Fumiaki Kitahara; writing—original draft preparation: Keiko Fukumoto; writing—review and editing: Tomohiro Nishizono, Fumiaki Kitahara; funding acquisition: Keiko Fukumoto; Resources: Keiko Fukumoto; supervision: Keiko Fukumoto.

Corresponding author

Correspondence to Keiko Fukumoto.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Handling editor: Ignacio J. Diaz-Maroto.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

Reprints and permissions