Contribution of tree community structure to forest productivity across a thermal gradient in eastern Asia

Synthetic data for Fig. 1

To provide examples of the proposed two hypotheses, i.e., species-response hypothesis and community structure hypothesis, for Fig. 1, we generated synthetic data assuming bivariate lognormal distributions of species relative woody productivity p_i and species standing biomass B_i, where i for species identity, with log-log linear, (or power-law) correlations, ln p_i = k + b ln B_i, as in left-hand panels of Fig. 1. The slope (scaling exponent) b is common at –0.15, and the constant k = –3.4 and –3.8 for tropical and temperate forests respectively for species response hypothesis (Fig. 1a), whereas k = –3.6 for both ‘tropical’ and ‘temperate’ forests for the community structure hypothesis (Fig. 1b). Mean ln B_i are –0.6 for two forests in Fig. 1a, while they are –1.0 and –0.2 for tropical and temperate forest respectively in Fig. 1b, Standard deviations of ln B_i and ln p_i are 2.0 and 0.65 respectively for all forests, except those in tropical forest in Fig. 1b are 1.6 and 0.6, respectively. In the left-hand panels, the B_i axis ranges 0.005–500 (Mg C ha^–1), and the p_i axis ranges 0.001–0.5 (yr^–1). In the right-hand panels, the axis for B = Σ_i B_i ranges 50–500 (Mg C ha^–1) and the axis for P = Σ_i p_i B_i ranges 0.5–20 (Mg C ha^–1 yr^–1).

Forest plot data

We selected 60 forest plots located in old-growth forests along the climatic gradient of insular eastern Asia, located on Java (3 plots), Kalimantan (5 plots), Peninsular Malaysia (2 plots), Taiwan (6 plots), and the Japanese archipelago (44 plots), ranging from 6.8°S to 44.4°N latitude and from 20 to 1,880 m in elevation (Supplementary Fig. 1, Supplementary Data 1). We collected climate data for all the plots for the period 1981–2010 from CHELSA version 2.1³⁹; these are the period-average annual and monthly ground surface mean temperature, precipitation, and potential evapotranspiration. The potential evapotranspiration was estimated by Hargreaves-Samani equation⁴⁰ based on monthly data of these climatic variables. Supplementary Data 2 presents mean annual temperature (MAT, °C), annual precipitation (AP, mm yr^–1), annual potential evapotranspiration (PET, mm yr^–1), monthly-data-based Thornthwaite moisture index (TMI) and the climatic types defined by TMI²⁶. The target region is in Asian monsoon climate^41,42, and moist forest ecosystems predominate from tropics in Southeast Asia to sub-boreal biomes in northern Japan. Across 60 plots, MAT ranges from 2.0 °C to 26.6 °C, AP-PET ranges from 58.5 to 5049 mm yr^–1, and plots are classified as “perhumid” or “humid” by TMI (Supplementary Data 2); the smallest TMI for the plot in cloud forest on Hahajima Island, oceanic Ogasawara Islands, where AP-PET was +217 mm yr^–1 (against +58.5 by CHELSA³⁹) based on the weather station records on the island. AP-PET sowed no correlation with MAT or with any forest structural or dynamic variable, in contrast to MAT exhibiting significant correlations to all forest variables (Supplementary Fig. 5). We therefore mainly employ MAT to quantify climatic dependence of the 60 plots. According to bioclimatic classification of the region^43,44, we define forest biomes into tropical (MAT ≥ 24 °C), subtropical (20–24 °C), warm-temperate (12–20 °C), cool-temperate (5–12 °C) and sub-boreal or subalpine (<5 °C). Dominant tree life forms are evergreen broad-leaved in tropical, subtropical, and warm-temperate forests, deciduous broad-leaved in cool-temperate, and evergreen coniferous in sub-boreal forests (Supplementary Data 3). We employed biomes to approximate forests into groups, while we used MAT values of each plot to analyse temperature dependence of forest stands and species populations.

As in Supplementary Data 1, 60 plots are ~1 ha in horizontal area (56 plots are exactly 1 ha, and others are 0.95, 0.98, 1.04, and 1.05 ha). Data include 41 (out of 44) plots across Japanese archipelago (code a in Supplementary Data 1; the original data is available at https://www.biodic.go.jp/moni1000/findings/data/)⁴⁵. Two cool-temperate mixed broadleaf-conifer forest plots (code c)²⁰, subalpine spruce-fir forest (code h, raw data originally uploaded here by T. S. Kohyama), two lowland mixed dipterocarp forest plots in West Kalimantan (code j)⁴⁶, three lowland heath or peat-swamp forests in Central Kalimantan (code f)^47,48,49, three tropical montane forest plots in West Java (code d)⁵⁰, two lowland mixed dipterocarp forest plots in Peninsular Malaysia (code i)^51,52, two subtropical karst forest plots (codes c)⁵³, two subtropical montane forest plots (code d)⁵⁴, and two warm-temperate montane forests in Taiwan^55,56. Where we have large continuous forest plots (codes b, c, g, i, h), we sampled every two 1-ha plots most distant from the other.

In each plot, tree censuses of stem diameter at breast height at 1.3 m above ground, or a marked position above buttress, were carried out for all stems ≥5 cm diameter appeared in either of two censuses with interval of ~5 years (ranging 2.5 and 8.5 years). We estimated individual-tree total aboveground carbon mass and leaf carbon mass from stem diameter using local allometric equations if available; otherwise, we used the generic allometric equations for eastern Asia⁵⁷ with species-specific wood densities^58,59. Generic estimates capture 89–99% of site-specific measures of tree mass⁵⁷. We used genus or family averages of wood density when species-specific data were not available, and when even these were not available, we used lifeform-specific equations (evergreen broadleaf, deciduous broadleaf, and evergreen conifer). We used the factor of 0.5 to convert oven-dry mass to carbon mass. Local equations are obtained in mixed dipterocarp forest in Peninsular Malaysia close to code i plots⁵¹ applied for these plots. The other from mixed dipterocarp forest in East Kalimantan⁶⁰ were applied to the same forest-type plot in West Kalimantan (code j) and montane forests in West Java (codes d) as well. Equations for heath forest and peat-swamp forest in Central Kalimantan⁴⁹ were applied to the same-site plots (code f) with respective forest types (‘Lahei-2’ as peat-swamp; other two as heath forest). In these local equations, we estimated tree height from stem diameter based on plot-specific inventory data, and obtained plot-specific coefficients of the extended allometry with asymptotic height that link diameter and height^49,50,51,61.

Species and forest productivity

We applied the following procedures²⁰ for estimating per-area forest aboveground woody productivity, i.e., aboveground net productivity by tree growth of surviving stems and ingrowth by recruited stems, that reduce the influence of inter-census interval and among-population variation²⁰. We estimate per-area aboveground biomass of species i in a plot, B_i (Mg C ha^–1), and species relative (i.e., per B_i) aboveground woody productivity p_i (yr^–1) as follows. For a census interval of T (yr), we obtained the aboveground biomass of a species i at the first (B_0i) and second census (B_Ti) as the sum of individual tree mass that were alive at each census, divided by the plot area (Mg C ha^–1). We also obtained the survived fraction of initial biomass (Bsurv_0i) as the sum of alive tree biomass at the first census that survived until the second census, divided by plot area. Then, using our methods²⁰, the estimated instantaneous relative aboveground woody productivity of species i is

and the period-mean aboveground biomass of species i over the two censuses is

The species absolute aboveground woody productivity is P_i = p_iB_i (Mg C ha^–1 yr^–1). We also obtain instantaneous relative aboveground woody loss rate (due to tree mortality) of species i by l_i = ln(B_0i/Bsurv_0i)/T, and absolute loss rate of i by L_i = l_i B_i, that counterbalance p_i and P_i, respecticvely²⁰. The forest period-mean aboveground biomass B (Mg C ha^–1) and aboveground woody productivity P (Mg C ha^–1  yr^–1) are respectively²⁰

and

Similarly, forest-level rate of aboveground woody loss is L = Σ_i L_i. In the provided code in the Zenodo repository (https://doi.org/10.5281/zenodo.7668416), we adopted our generalised estimation scheme for relative woody productivity p and relative woody loss rate l by tree mortality of populations for varied inter-census intervals among individual trees²¹. We define forest tree abundance N_i (ha^–1), or per-area stem count (≥5 cm stem diameter) of population i to be the sum of population period-mean per-area stem counts⁶¹.

In each plot, we selected every species (or morpho-type) with two or more trees that survived through two censuses, and combined all other rarer species (<two surviving trees per plot) within an aggregated multi-species population. Total number of species across 60 plots was 1587, that of per-plot species populations was 3807, and that without aggregated rare populations was 2604. When any species ≥ two surviving trees showed non-positive p_i, due to diameter decrease of large tree(s) during the corresponding period, we discarded those populations as in p_i–B_i model fitting in Fig. 2a, b (53 out of 2604 populations, i.e., 2.1%), but we included all species with aggregated rare species and with p_i ≤ 0 in forest-level woody productivity and loss rate estimates.

Productivity model fitting

To fit and quantify power-law models, we applied linear models with log-transformed variables assuming that response variables are lognormally distributed such as

where j is plot (or species population) identity, y_j is a response variable, x_j and MAT_j are explanatory variables, μ_j and σ are j-specific mean and common residual deviation, and k, b and c are model coefficients to be estimated. The logarithmic model of Eq. (6) is converted to the power-law model,

where a = exp(k + σ²/2) by adjusting the mean of lognormally distributed y_j from normally distributed ln y_j. In species population level model fitting, we excluded ‘aggregated’ rare species and records with non-positive p_i as y_j, but we included these in forest plot-level analyses. We present conditional R² values⁶² for log-log linear mixed models of species aboveground biomass against species maximum tree size or species abundance (per-area tree count) with plot-specific constant terms in Supplementary Fig. 2. Chemical theory has stimulated that MAT-dependence of turnover rates is related to the inverse of absolute temperature (°C), i.e., 1/(273.15 + MAT)^22,30. We used a simpler formulation of Eq. (7) over our limited range of MAT values (0–30 °C), noting the near linear dependence between MAT and 1/(273.15 + MAT) (R² = 0.999).

We carried out data analyses and graphical presentations with R version 4.0.5⁶³ and Python version 3.9⁶⁴. All data at the level of individual trees and species populations, and the R code used to generate Fig. 4 are provided in the Zenodo repository (https://doi.org/10.5281/zenodo.7668416).

Null model data of temperature dependence

To test the species response hypothesis (Fig. 1b) separated from community structure difference among plots, we calculated woody productivity of each species population responding to the plot location in terms of the plot-specific power-law constant, a_plot, of species productivity (p_i) versus biomass (B_i), i.e., p_i ~ a_plot B_i^b (cf. Figure 2a), and that frequency distribution of species biomass is the same as the distribution of all B_i’s across all the 60 plots. Our species response estimate of forest-level woody productivity for each plot is:

This estimate assumes the plot species richness is proportional to the plot biomass B. In this estimation procedure, we left the aggregated population of rare species unchanged.

To disentangle the contribution of species response hypothesis and community structure hypothesis (Fig. 1), we generated two null-model woody productivity data of 60 plots. To represent the effect of community structure on woody productivity (Fig. 1a), we generated “replaced data” from the observation data, in which we replaced the productivity of a given species i, P_i, with the value of different species populations with similar species biomass B_i. To do this, we performed a weighted random sampling 10,000 times from the species population pool (excluding the target species) of all plots, using the inverse of the square of the difference in ln B_i between species populations as the weight. We took the mean value of the resampled data as the “replaced” P_i for each species population i. The aggregated population of rare species were left unchanged, and excluded from this replacement procedure. The sum of these values was taken as plot-level productivity estimate representing the community structure hypothesis, P_commStr, in Fig. 4b. Supplementary Fig. 6 compares species woody productivity between original and replaced data.

To quantify the contribution of our two hypotheses, we applied regression model of ln P_specRes (Fig. 4b) and that of ln P_commStr (Fig. 4c) to the original data with respect to ln P, and compared their residual variances to that of the original data model (Fig. 4a) as: [variance minus the residual variance from projected model estimates] divided by [variance minus the residual variance from original model estimates]. This ratio is 1 if a null model completely explains ln P as original model does, and is 0 if the null model explains nothing of the variance with respect to ln P.

Quantifying biases of standard productivity estimates

Standard, ‘simple’ estimates of aboveground woody productivity, P_simple (Mg C ha^–1 yr^–1), is the sum of absolute aboveground mass gain by survived and recruited trees of any species, i.e., P_simple = Σ_i(B_Ti−Bsurv_0i)/T, setting the initial mass for recruited trees is the mass at the threshold tree size (5 cm stem diameter here)¹⁷. We have theoretically shown that this conventional woody productivity estimate of a population i, P_simple_i is influenced by inter-census interval T, latent relative woody productivity p_i (Eq. (1)) and loss rate l_i, thus by inter-specific variation of p_i and l_i, in such manner that P_simple_i of species i with larger p_iT and (p_i−l_i)T are more largely underestimated²⁰. It is also influenced by the treatment of unrecorded initial tree mass^19,20.

We obtained standard, simple estimates of forest woody productivity of the 60 plots to evaluate them in comparison to our instantaneous estimates (Eqs. (1)–(4))²⁰. The simple productivity estimates were in average 85% in terms of species population woody productivity and 90% of our forest woody productivity, with large inter-specific and among plot variations (Supplementary Fig. 7a, b). At the species-level, simple estimates predicted the exponent of relative productivity versus biomass power-law at –0.077 ± 0.009 in contrast to that of –0.14 ± 0.01 in our estimates (Supplementary Fig. 7c, compared to Fig. 2a). Forest-level model of P_simple dependence on standing biomass and temperature was similar to our model of P, besides that P_simple was smaller than P by about 10% (Supplementary Fig. 7d, compared to Fig. 4a).

Aboveground net primary productivity

There were monthly records of litter fall collected using litter traps during the corresponding period in 22 out of our 60 plots. We estimated aboveground net primary productivity (i.e., the estimated aboveground woody productivity P plus recorded monthly-sum fine litter fall (i.e., canopy productivity⁹) during the same census period) to validate our analysis for this more inclusive measure of aboveground net productivity (ignoring non-structural organic carbon productivity). Our estimates of forest-level aboveground woody productivity P are log-log linearly correlated with aboveground net primary productivity (R² = 0.88; Supplementary Fig. 8a). Our model of forest woody productivity as functions of forest biomass and mean annual temperature (Fig. 4a) also apply to forest net primary productivity (Supplementary Fig. 8b).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Source: Ecology - nature.com