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 pi and species standing biomass Bi, where i for species identity, with log-log linear, (or power-law) correlations, ln pi = k + b ln Bi, 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 Bi 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 Bi and ln pi 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 Bi axis ranges 0.005–500 (Mg C ha–1), and the pi axis ranges 0.001–0.5 (yr–1). In the right-hand panels, the axis for B = ΣBi ranges 50–500 (Mg C ha–1) and the axis for P = ΣpBi 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.139; these are the period-average annual and monthly ground surface mean temperature, precipitation, and potential evapotranspiration. The potential evapotranspiration was estimated by Hargreaves-Samani equation40 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 TMI26. The target region is in Asian monsoon climate41,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 CHELSA39) 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 region43,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 Two cool-temperate mixed broadleaf-conifer forest plots (code c)20, 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)46, 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)50, two lowland mixed dipterocarp forest plots in Peninsular Malaysia (code i)51,52, two subtropical karst forest plots (codes c)53, two subtropical montane forest plots (code d)54, and two warm-temperate montane forests in Taiwan55,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 Asia57 with species-specific wood densities58,59. Generic estimates capture 89–99% of site-specific measures of tree mass57. 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 plots51 applied for these plots. The other from mixed dipterocarp forest in East Kalimantan60 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 Kalimantan49 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 height49,50,51,61.

Species and forest productivity

We applied the following procedures20 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 variation20. We estimate per-area aboveground biomass of species i in a plot, Bi (Mg C ha–1), and species relative (i.e., per Bi) aboveground woody productivity pi (yr–1) as follows. For a census interval of T (yr), we obtained the aboveground biomass of a species i at the first (B0i) and second census (BTi) 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 (Bsurv0i) as the sum of alive tree biomass at the first census that survived until the second census, divided by plot area. Then, using our methods20, 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 Pi = piBi (Mg C ha–1 yr–1). We also obtain instantaneous relative aboveground woody loss rate (due to tree mortality) of species i by li = ln(B0i/Bsurv0i)/T, and absolute loss rate of i by Li = li Bi, that counterbalance pi and Pi, respecticvely20. The forest period-mean aboveground biomass B (Mg C ha–1) and aboveground woody productivity P (Mg C ha–1  yr–1) are respectively20






Similarly, forest-level rate of aboveground woody loss is L = Σi Li. In the provided code in the Zenodo repository (, 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 trees21. We define forest tree abundance Ni (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 counts61.

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 pi, due to diameter decrease of large tree(s) during the corresponding period, we discarded those populations as in piBi 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 pi ≤ 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

$${{{{{rm{ln}}}}}},{y}_{j} sim {{{{{rm{Normal}}}}}}({mu }_{j},,sigma ),$$




where j is plot (or species population) identity, yj is a response variable, xj and MATj 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,

$${y}_{j}=a,{x}_{j}^{b}{{exp }}(c,{{{{{{rm{MAT}}}}}}}_{j}),$$


where a = exp(k + σ2/2) by adjusting the mean of lognormally distributed yj from normally distributed ln yj. In species population level model fitting, we excluded ‘aggregated’ rare species and records with non-positive pi as yj, but we included these in forest plot-level analyses. We present conditional R2 values62 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) (R2 = 0.999).

We carried out data analyses and graphical presentations with R version 4.0.563 and Python version 3.964. 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 (

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, aplot, of species productivity (pi) versus biomass (Bi), i.e., pi ~ aplot Bib (cf. Figure 2a), and that frequency distribution of species biomass is the same as the distribution of all Bi’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, Pi, with the value of different species populations with similar species biomass Bi. 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 Bi between species populations as the weight. We took the mean value of the resampled data as the “replaced” Pi 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, PcommStr, 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 PspecRes (Fig. 4b) and that of ln PcommStr (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(BTiBsurv0i)/T, setting the initial mass for recruited trees is the mass at the threshold tree size (5 cm stem diameter here)17. We have theoretically shown that this conventional woody productivity estimate of a population i, P_simplei is influenced by inter-census interval T, latent relative woody productivity pi (Eq. (1)) and loss rate li, thus by inter-specific variation of pi and li, in such manner that P_simplei of species i with larger piT and (pili)T are more largely underestimated20. It is also influenced by the treatment of unrecorded initial tree mass19,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))20. 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 productivity9) 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 (R2 = 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 -

Larval rockfish growth and survival in response to anomalous ocean conditions

When legislation to protect wildlife becomes a problem