The relative importance of abiotic and biotic factors in shaping forest biomass stocks and fluxes remains a controversial issue. Here, using data gathered from 39 1-ha plots located in flooded and terra firme mature tropical lowland forests of the Amazon and Orinoquia regions of Colombia, we evaluated the importance of climate, soil fertility, flooding as well as tree taxonomic/phylogenetic diversity and forest structural properties, in determining the aboveground biomass stocks (AGB; Mg ha-1) and aboveground woody productivity (AWP; Mg ha-1 y-1). Using information-theoretic multimodel inference and variance partitioning we found that forest structural features, such as the number of trees with diameter at breast height ≥ 70 cm and wood density, are the main drivers of variation in AGB. However, taxonomic diversity also contributes to AGB because it is associated with more large trees in these forests. In contrast, the key drivers of AWP in these forests were soil P and Mg concentration, with no significant effects of diversity indices. These findings emphasize the need to include other major soil cations than N and P (e.g., Mg) in experimental studies to improve our understanding about the extent to which soil fertility can modulate increases in forest AWP due to climate change. Terra firme forests had higher AGB stocks than flooded forests, but both had similar AWP; and we found similar results for the drivers of AGB and AWP between flooded and terra firme forests. Our results provide limited evidence for strong effects of plant diversity on AGB or AWP. Therefore, we call for caution on generalizations of nature-based initiatives aiming to preserve diversity based on maximizing carbon stocks and productivity due to the complex nature of the processes controlling carbon accumulation and carbon fluxes in tropical forests.

This study was conducted using tree census data collected at 39 1-ha forest inventory plots situated across a latitudinal range from -3.8° to 4.8°, a longitudinal range from -77.1° to -67.8°, and an elevation range from 79 m asl to 963 m asl, in Colombia.  The mean annual temperature (MAT) of plots ranged from 21.5 to 25.5 °C (mean = 24.7 ± 0.6 °C; mean ± SD) and mean annual precipitation (MAP) of the plots ranged from 1797 to 3716 mm y^-1 (3076.4 ± 299.9 mm y^-1) (External Databases S1). There were 8 plots located in the Orinoquia region (González-Abella et al. 2021) and 31 in the Amazonia region of Colombia (Fig. S1). Plots were censused between 2005 and 2021. In the locality of Amacayacu, the data employed come from a continuous 25-ha plot located on Tierra firme soils derived from the Pebas formation (Zuleta et al. 2017).  Therefore, to avoid pseudo replication and over representation of this site, we only used the minimum, maximum, and mean value of the AGB, AWP, and AWRT at the 1-ha scale. Including these three values allow us to have a landscape representation of the variation of carbon stocks and fluxes in this geological formation characterized for having the relatively more fertile soils in the Amazon Tierra firme forests (Hoorn 1994). In total, 11 out of the 39 plots were located on floodplains, 8 of which were in Igapo and 3 in Varzea.

In each plot, all stems of trees and palms (hereafter trees) with tree diameter at breast height (DBH; tree diameter at 1.3 m height) ≥ 10 cm were measured. We corrected all DBH measured over 1.4 m height using the trunk taper equation proposed by Cushman et al. (2021). Some trees with an annual negative growth in DBH < 1 cm, were corrected to a negative annual growth of DBH equals to - 1.0 cm. 

The aboveground biomass (AGB) of each tree was estimated using the allometric equation proposed by Chave et al. (2014), defined as: ; where AGB (kg) is the estimated aboveground biomass, DBH (cm) as defined above, H (m) is the estimated total height, and WD (g cm^-3) is the stem wood density (see Supplementary Information for more details). The AGB of each palm was estimated using the allometric equation proposed by Goodman et al. (2013) defined as:  where dmf is the dry mass fraction calculated as the ratio between dry mass and fresh mass. AGB per plot was estimated as the sum of AGB of all individuals in each plot.

All plots were censused at least twice (elapsed time ranged between 2 and 10 years), and the aboveground woody productivity (AWP in Mg ha^-1 y^-1), and aboveground woody residence time (AWRT in y) of each plot were estimated. The AWP was calculated as the sum of the annualized AGB recruitment and AGB growth. The annualized recruitment AGB (rAGB; Mg ha^-1 y^-1) is the sum of the AGB of individuals recruited into DBH ≥ 10 cm between censuses divided by the time between measurements. However, for each tree recruited (DBH ≥ 10 cm), we subtracted the corresponding AGB associated with a tree of 9.99 cm (i.e., just below the detection limit) to avoid overestimations of the overall increase in AGB due to recruitment (Talbot et al. 2014). The annualized AGB growth (gAGB; Mg ha^-1 y^-1) was estimated as the sum of the increase in AGB of all individuals with DBH ≥ 10 cm that survived between censuses divided by the time between censuses. The AWRT was calculated as the quotient between the annualized values of AGB and AWP in each plot (Muller-Landau et al. 2021). We also calculated the net change in AGB as the annualized difference in AGB between the final and first census of each plot divided by the elapsed time between measurements (Duque et al. 2021). 

The aboveground biomass (AGB) of each tree was estimated using the allometric equation proposed by Chave et al. (2014), defined as: ; where AGB (kg) is the estimated aboveground biomass, DBH (cm) as defined above, H (m) is the estimated total height, and WD (g cm^-3) is the stem wood density (see Supplementary Information for more details). The AGB of each palm was estimated using the allometric equation proposed by Goodman et al. (2013) defined as:  where dmf is the dry mass fraction calculated as the ratio between dry mass and fresh mass. AGB per plot was estimated as the sum of AGB of all individuals in each plot.

All plots were censused at least twice (elapsed time ranged between 2 and 10 years), and the aboveground woody productivity (AWP in Mg ha^-1 y^-1), and aboveground woody residence time (AWRT in y) of each plot were estimated. The AWP was calculated as the sum of the annualized AGB recruitment and AGB growth. The annualized recruitment AGB (rAGB; Mg ha^-1 y^-1) is the sum of the AGB of individuals recruited into DBH ≥ 10 cm between censuses divided by the time between measurements. However, for each tree recruited (DBH ≥ 10 cm), we subtracted the corresponding AGB associated with a tree of 9.99 cm (i.e., just below the detection limit) to avoid overestimations of the overall increase in AGB due to recruitment (Talbot et al. 2014). The annualized AGB growth (gAGB; Mg ha^-1 y^-1) was estimated as the sum of the increase in AGB of all individuals with DBH ≥ 10 cm that survived between censuses divided by the time between censuses. The AWRT was calculated as the quotient between the annualized values of AGB and AWP in each plot (Muller-Landau et al. 2021). We also calculated the net change in AGB as the annualized difference in AGB between the final and first census of each plot divided by the elapsed time between measurements (Duque et al. 2021)

Climate variables at each plot location were extracted from the CHELSA (Karger et al. 2017) bioclimatic rasters at a resolution of 30-arcsec (approximately 1 km² at the equator). We first extracted the 11 climatic variables associated with temperature and the 8 climatic variables associated with precipitation available in the CHELSA repository (see Supplementary Information). Then, we selected mean annual Temperature (MAT; °C), total annual precipitation (TAP; mm y^-1), and precipitation seasonality (PS; kg m^-2), defined as the coefficient of variation of the monthly precipitation estimates expressed as a percentage of the mean of those estimates (i.e. the annual mean). Once these three variables were selected, the remaining 16 variables included showed a correlation ≥ 0.7 with MAT, TAP, or PS and were excluded.

To estimate soil fertility, samples of soil A horizon (i.e., the mineral soil after removing the organic layer) were collected from a minimum of five points in each plot at a depth of 10–30 cm. The five samples from each plot were then combined, and a 500g composite sample was taken and air-dried after removing macroscopic organic matter. Sodium (Na; mg kg^-1), aluminum (Al), Cation exchange capacity (CEC), total phosphorus (P; mg kg^-1), calcium (Ca; mg kg^-1), magnesium (Mg; mg kg^-1), potassium (K; mg kg^-1), and contents of organic matter (OM %) were analyzed for each plot (see Supplementary Information for details about lab analysis).

Before assessing phylogenetic and taxonomic diversity, all species and genus names were checked and standardized using the Taxonomic Name Resolution Service (Boyle et al. 2013). In the dataset, 90.4% of stems were identified to species level, 5.3% to genus, 1.0% to family, and 3.3% remained unidentified. Based on the individuals fully identified to the species level, we assessed the number of species in each plot as a metric of taxonomic richness (see Duque et al. 2017), and the inverse of the Simpson index as a metric of taxonomic diversity (see Coelho de Souza et al. 2019). These metrics were estimated using vegan R package (Oksanen et al. 2022).

We calculated three metrics of phylogenetic diversity: phylogenetic diversity sensu stricto (PD), Net Relatedness Index (NRI) and the Nearest Taxon Index (NTI). NRI and NTI were weighted by abundance. To calculate these indexes, we generated a phylogenetic tree from all fully identified species using the backbone phylogeny proposed by Smith and Brown (2018); this analysis was done employing V.PhyloMaker (Jin and Qian 2019). The PD of each plot was calculated as the total sum of the phylogenetic branch lengths connecting the co-occurring species in each plot along the minimum spanning path to the root of the tree. The NRI and NTI are based on the mean pairwise distance and the mean nearest pairwise distance, respectively (Webb 2000). These observed distances were standardized with respect to random distribution obtained from a null model. We used an independent swap null model that allow us to fit species frequency and abundance in each randomization. The standardized values are multiplied by -1. NRI and NTI positive values represent phylogenetic clustering and negative values represent phylogenetic overdispersion. These metrics were estimated using the picante package in R (Kembel et al. 2010).

Forest structural properties were represented by forest stand descriptors, such as number of individuals (N), average wood density per plot (WD), number of trees with DBH ≥ 70 cm (D70), and maximum DBH (Dmax), which have shown to play an important role in shaping AGB, AWP, and AWRT in tropical forests (Slik et al. 2013, Fauset et al. 2015, Duque et al. 2021). Other biological processes such as biological associations (e.g., micorrhizas) and species competition are not considered here.

We used an Information Theory (IT) natural model-averaging technique (Burnham and Anderson 2002) to identify the dominant ecological drivers of AGB, AWP, and AWRT across the Amazonia and Orinoquia regions of Colombia. The explanatory variables associated with our hypotheses were: climate and soil fertility as abiotic factors, and phylogenetic/taxonomic diversity and forest structure as biotic factors. The natural model-averaging technique allows us to infer and predict AGB, AWP, and AWRT based on a set of alternative models. The natural model-averaging technique calculates the average of each variable´s parameter estimates over the models where the variable was selected (Symonds and Moussalli 2011, Galipaud et al. 2017). The explanatory variables were previously standardized to have mean = 0 and standard deviation = 1, and then each parameter standardized by their partial standard deviations (Cade 2015). The use of the partial standard deviations (Cade 2015) aims to correct for the likely effect of multiple correlation among variables. The partial standard deviation (s*ₓᵢ) is calculated as follows.