Tree aboveground biomass and species richness of the mature tropical forests of Darien, Panama, and their role in global climate change mitigation and biodiversity conservation

The remote forests of the Darien region in eastern Panama are among the last remnants of relatively undisturbed forest habitat in the Central American isthmus. Despite decades of efforts by the government, nongovernmental organizations, and civil society, including Indigenous peoples, to protect the region's natural heritage, it remains under significant threat due to widespread illegal logging. Now, the Panamanian government is considering the mechanism, Reducing Emissions from Deforestation and Forest Degradation (REDD+), as another option to limit forest loss. Central to the proper functioning of REDD+ is the need to reduce uncertainties in estimates of aboveground biomass (AGB). These estimates are used to establish realistic reference levels against which additional contributions to reducing carbon dioxide emissions from the loss and degradation of forests can be financially compensated. Also, highly desirable to REDD+ is the achievement of biodiversity cobenefits. REDD+ investments will likely be directed primarily to areas where the potential to simultaneously mitigate climate change and conserve biodiversity is highest. Here, we present the results of a field‐based forest carbon inventorying method tested in Darien's mature forests with the participation of Embera and Wounaan Indigenous peoples. We also explore whether variations in field‐based estimates of AGB across mature forests, in both undisturbed and disturbed areas, are detectable through free and readily‐available remote sensing data sources. Furthermore, we examine and compare AGB and tree species richness in Darien with other well‐studied forest sites across the tropics. Our findings reveal that Darien's forests play a crucial role globally and regionally in storing carbon and housing biodiversity, and support the imperative need to protect these forests in a culturally appropriate manner with the region's Indigenous peoples.


| INTRODUCTION
and Ramsar Wetland of International Importance. Most recently, the Government of Panama identified Reducing Emissions from Deforestation and Forest Degradation (REDD+) as an additional promising option to stem the loss of forests in Darien, and elsewhere in Panama (Melgarejo, Corro, Ruiz-Jaén, Calderón, & Sánchez de Stapf, 2015). Since 2007 (United Nations Framework Convention on Climate Change, 2008), REDD+ has been at the forefront of international policy efforts to curb atmospheric carbon dioxide (CO 2 ) emissions from land use and land-use change. Given that deforestation and forest degradation account for 9% of CO 2 emissions globally (Le Quéré et al., 2015), REDD+ could play a preponderant and cost-effective role (Stern, 2007) in mitigating climate change.
The present paper is part of a participatory project that developed a culturally appropriate method to empower Embera and Wounaan Indigenous peoples, the original custodians of Darien's forests, to estimate forest carbon stocks (Mateo-Vega et al., 2017). This project was requested by Indigenous authorities aware of a previous study by Asner et al. (2013), which pooled most of Darien's mature forests into a single, very high carbon density class of ≥130 Mg/ha (i.e., 277 Mg/ha of aboveground biomass-AGB) using airborne light detection and ranging (LiDAR) technology. Results contrasted widely with those from Panama's field-based National Forest and Carbon Inventory F I G U R E 1 Location of undisturbed (green) and disturbed (red) mature forest plots in Darien (Melgarejo et al., 2015), which found that mature forests, such as those of Darien, only contained~80 Mg ha of carbon (i.e., 170 Mg/ha of AGB). Indigenous authorities were interested in quantifying forests carbon stocks using fieldbased measurements to validate the REDD+ potential of their forests, and engage in informed discussions with REDD+ proponents.
Accordingly, here we analyze sources of AGB variation in 30 plots of approximately 1 ha each, distributed across a large mature forest landscape. We hypothesize that (a) the relative density of large versus small trees, (b) forest types, and (c) disturbance resultant from traditional forest use could explain heterogeneity at the landscape level. As field inventories are expensive and time consuming, we also assess whether spatially coarse remote sensing-derived products capture differences in AGB between undisturbed and disturbed mature forest areas based on forest cover (Sexton et al., 2013) and tree height (Simard, Pinto, & Fisher, 2011). To fully capture the conservation value of Darien's forests, we then compare our AGB estimates and tree species richness with other well-studied forests across the tropics.

| Study site and data collection
Previous studies (Mascaro et al., 2011) have questioned the statistical power of single large plots in capturing landscape-level AGB. Therefore, as reported in Mateo-Vega et al. (2017), we chose to establish 30, square-shaped, nested plots of approximately 1 ha, with four internal 12 × 12-m subplots, in forests within five distant Indigenous territories that are part of Tierras Colectivas Embera y Wounaan of Darien (Figures 1 and 2). The realized average plot size, given topography and human error, was 0.93 ha, and the total area sampled was~28 ha (Mateo-Vega et al., 2017). The four subplots within each plot correctly measured 144 m 2 .
We distributed our plots in clusters of 5-8 across 3,500 km 2 (~5% of Panama's territory) (Figure 1), establishing 18 in Tropical Moist Forest, 8 in Pre Montane Wet Forest, and 4 in Tropical Wet Forest, following the Holdridge Life Zone classification system (Holdridge, 1967) ( Figure 2). Twenty-five plots were established in territories that are not accessible by road, and that require 6-20 hr to reach via boat, dugout canoe and/or hiking, after a 7-hr drive to Puerto Quimba, entry point to the easternmost part of Darien. The five remaining plots were located in a territory accessible by road, but require a 2-5 hr hike to reach them, after a 6-hr drive from Panama City. The sites and number of plots established per territory was determined in consultation with the Indigenous authorities. Fifteen plots were established in mature forests with evidence of disturbance (e.g., tree stumps, trails, and nearby agriculture; Mateo-Vega et al., 2017) resultant from traditional land uses (e.g., selective logging, harvest of nontimber forest products), and the remainder in undisturbed areas (Figures 1 and 2 of sampling design, we recognize subplots to be nested within plots, and plots to be nested under territories, as well as under forest types. In each territory, plots were established in intact and disturbed forests, and therefore these factors can be adequately considered crossed. The diameter-at-breast-height (DBH) and height of all living large trees ≥50 cm DBH (including palms and excluding lianas) were measured in each 1-ha plot, while all small trees ≥10 cm to <50 cm DBH were measured in four internal subplots (Mateo-Vega et al., 2017). We measured the DBH and height of 1,401 living trees, and identified these in the field or from voucher specimens at the species (144), genus (58), or family (64) level. All but 72 trees were identified at one of these levels. Because tree wood density is a key predictor of AGB , the wood density for each tree was determined using the Global Wood Density database Zanne et al., 2009). For trees that professional botanists were unable to identify at any level, a regional wood density value for Central America was assigned Zanne et al., 2009).

| Tree biomass, canopy height, and species richness variation
As reported in Mateo-Vega et al. (2017), AGB for each tree was estimated using the Chave et al. (2014) pantropical allometric model. For palms, we used the Goodman et al. (2013) allometric model. We singled out three sources of variation in AGB hypothesizing that the relative density of large versus small trees could explain heterogeneity at the landscape level, and that large trees density could result from differences across forest types, natural forest dynamics, or human intervention.
Despite representing a small percentage of trees in a forest, large trees account for a sizeable proportion of total AGB and may be key predictors of stand-level biomass (Slik et al., 2013). To elucidate sources of large tree AGB variation, and recognizing incomplete nesting of factors ( Figure 2), we ran separate one-way analysis of variance (ANOVA) tests to compare the plot-level mean AGB of large trees among the five territories, three forest types, and three forest types while controlling for undisturbed and disturbed areas. Due to imbalanced sampling, we ran bootstrapping tests (10,000 iterations) to select the minimum number of plots for each category. We also compared plotlevel mean AGB of large trees between undisturbed and disturbed plots using a Welsh two-sample t-test. We used AGB values scaled to 1 ha based on plot measurements for all tests.
To elucidate AGB variation for small trees within and among territories, we ran a two-level nested ANOVA on subplot mean AGB. The model employed subplots nested within plots in each territory (first level), and then subplots in plots nested among the five territories (second level). Due to imbalanced sampling, we ran a bootstrapping test (10,000 times) to randomly select five plots per territory and three subplots per plot (i.e., minimum number of plots and subplots sampled per territory and plot, respectively) for each iteration, thus restoring equal sample size.
We also examined AGB variation at the landscape level by scaling mean small-tree AGB in each plot to a 1-ha area, then summing these values with large tree AGB to obtain total AGB per ha. An ANOVA was then carried out to test for significant differences in plot-level AGB among the five territories using a bootstrapping test (10,000 iterations) to randomly select five plots per territory. We examined differences between plot-level AGB in undisturbed and disturbed forest plots using a Welsh two-sample t-test. We further used the coefficient of variation (CV) of all plots, and among or within the five territories, three forest types, and two forest states, to quantify variation. All statistical analyses were conducted in R Version 3.3.1 (R Core Team, 2016).
Canopy structure characteristics are key for estimating AGB with field-based (Chave et al., 2005) and remotesensing methods (Asner et al., 2013). To test for significant differences between undisturbed and disturbed forest plots, we compared our mean large-tree height field measurements using a Welsh two-sample t-test. To corroborate if our upper canopy height measurements are consistent with remotelysensed data, we compared them with the Global Forest Canopy Height 2005 dataset (Simard et al., 2011). Due to cloud cover in the satellite images, we excluded one undisturbed and three disturbed forest plots from our analysis. Furthermore, we tested for significant differences in height measurements between undisturbed and disturbed forest plots based on the satellite images using a Wilcoxon/Kruskal-Wallis test. We also examined differences in tree cover, an indicator of forest degradation, for undisturbed and disturbed forest plots based on the freely available Landsat vegetation continuous field (VCF) satellite imagery layers of Sexton et al. (2013) for 2010 and 2015. We excluded eight and six plots for both time periods, respectively, due to cloud cover. These layers estimate the percentage of horizontal ground covered by woody vegetation greater than 5 m in height, at a 30 m spatial resolution (Supporting Information Data S1). We tested for significant differences in percentage tree cover for undisturbed and disturbed forest plots using a Wilcoxon/Kruskal-Wallis test.
Finally, to examine tree species variation across the five territories, we measured beta diversity (Sorensen), including species turnover (species replacement at one site by others) and nestedness (absence and nonreplacement of species from one site to another) (Baselga & Orme, 2012). Each territory, and plots contained therein, was considered a single site. We examined various beta diversity components for (a) all living trees, (b) large and small trees, (c) trees in undisturbed and disturbed forest plots, (d) trees in different forest types (excluding Tropical Wet Forest, found only in one territory), and (e) combinations of large and small trees with undisturbed and disturbed sites, and forest types (excluding all Tropical Wet Forest and Pre Montane Wet Forest plots in disturbed forest areas, found in only one territory). Bootstrapping tests of 1,000 iterations were carried out to randomly select the number of territories, and plots per territory, necessary to run the comparisons across all categories and permutations. We used the betapart function in R by Baselga and Orme (2012), namely the multi-sites dissimilarities computation, to carry out the analyses.

| Establishing the importance of Darien's AGB and tree diversity in relation to other tropical forests
To compare Darien to other well-studied forests throughout the Neotropics, tropical Asia and Africa, we used AGB estimates of our plots with the same allometric models and wood density values (Chave et al., 2005Zanne et al., 2009), as synthesized in Réjou-Méchain et al. (2014). The comparison only considered undisturbed sites, with the exception of the Luquillo Forest Dynamics Plot (Puerto Rico) (Thompson et al., 2004).
Comparing tree species richness from our plots with other sites is challenging given differences in plot size and configuration (i.e., single large plot versus numerous 1-ha plots across a landscape). Nonetheless, we compared our plots with Barro Colorado Island (BCI) in central Panama, a reference point for tree diversity in Panama and the Neotropics, more generally. For this, we ran a bootstrapped analysis (1,000 iterations) on the 2010 Census 7 BCI data (Condit, 1998;Hubbell et al., 1999;Hubbell, Condit, & Foster, 2010). The 50-ha plot was divided into 50, 1-ha plots with their respective subplots as described in Mateo-Vega et al. (2017). For each iteration, 30 plots were randomly selected, and the mean number of species estimated for large trees at the plot level, and small trees at the subplot level.

| AGB and tree species variation
We found no significant difference in mean per-ha AGB of large trees among territories (F 4,20 = 1.99, p value = 0.23) or among the three forest types sampled (F 2,12 = 0.91, p value = 0.60). Likewise, large tree AGB per ha among the three forest types in undisturbed forest plots, exhibited no significant difference (F 2,6 = 2.51, p value = 0.21). The same analysis for disturbed forest plots, only comparing Pre Montane Wet Forest and Tropical Moist Forest (i.e., there was only one Tropical Wet Forest plot in disturbed forests), did not reveal significant differences either (F 1,1 = 5.44, p value = 0.45). However, large tree AGB ranged from 121 to 349 Mg/ha in undisturbed forest plots, and from 12 to 200 Mg/ha in disturbed forests, and the difference between the two forest states was significant (t = −5.751, p value = 3.603e−06).
We also examined variation in small-tree AGB using a nested ANOVA within and among territories. Mean subplot small-tree AGB among the five territories ranged from 0.96 to 1.87 Mg/ha. The results of the ANOVA indicate no significant variation for subplots nested within plots in individual territories (F 20,50 = 0.79, p value = 0.70), nor among the five territories (F 4,20 = 3.3, p value = 0.07). The absence of variation confirmed that our scaling-up approach, applying small-tree biomass to the entire plot, was adequate.
For all 30 plots, mean large tree AGB accounts for 63.3% of total AGB, and is found in 30 (SD 14) trees/ha on average ( Figure 3). Small-tree AGB accounts for 36.7% and is found on average in 374 (SD 112) trees/ha. In undisturbed forest plots, the proportion of mean AGB for large trees increases to 68.5% (mean of 39 [SD 11] trees/ha) and lowers to 31.5% (mean of 369 [SD 128] trees/ha) for small trees. In disturbed forest plots, mean AGB of large trees decreases to 55.3% (mean of 22 [SD 11] trees/ha), thus increasing small-tree AGB to 44.7% (mean of 337 [SD 97] trees/ha).
Total AGB among all 30 plots ranged from 87 to 460 Mg/ha with a mean of 283 (SD 89) Mg/ha (Figure 4). This is equivalent to 133 Mg of carbon/ha, similar to the value assessed by Asner et al. (2013)-i.e., 130 Mg carbon/ha-for the forests of Darien using airborne LiDAR (2% difference). Mean AGB among the five territories ranged between 206 and 366 Mg carbon/ha, and we found no significant difference among these (F 1,2 = 5.3, p value = 0.05). In undisturbed forest areas, mean plot-level AGB was 54% greater than in disturbed forest areas. We found significant differences between plot-level AGB for undisturbed and disturbed forest areas (t = −4.984, p value = 2.931e −05) across the five territories. Total mean AGB for plots in Tropical Moist Forests in undisturbed areas was highest (379 [SD 51] Mg/ha), and also exhibited the largest difference (71%) with plots in disturbed forest areas. In the case of Pre Montane Wet and Tropical Wet Forest plots, the differences in mean AGB between undisturbed and disturbed areas was also large, that is, 47 and 32%, respectively.
The CV allowed us to quantify the amount of variation in AGB across the landscape. We found that ecological determinants of forest carbon stocks (e.g., forest types) cause less variation in forest AGB than anthropogenic activity, since the CV among forest states (i.e., undisturbed and disturbed) amounted to 30%, whereas only 7% amongst forest types (Table 1).
To test for significant differences in canopy height and tree cover between undisturbed and disturbed forest plots, we used field measurements and freely available satellite data. We focused on these two categories for the satellite imagery analysis to corroborate our field measurements, which yielded significant differences. The mean canopy height of large trees in disturbed forests based on field measurement was 30.3 m (SD 8) and 31.8 m (SD 8) for undisturbed forests. The results were statistically significant (t = −2.459, p value = 0.014). The analysis of upper canopy tree height with the satellite-based Global Forest Canopy Height dataset (Simard et al., 2011), revealed that undisturbed forest plots had a mean height of 33.2 m (SD 6), while disturbed forest plots had a mean height of 29.6 m (SD 5). The difference was not statistically significant. When comparing percentage total tree cover between disturbed and undisturbed forest plots using the satellite-based VCF forest cover layers of Sexton et al. (2013), no significant differences were found (Z = 0.11577, p value = 0.908).
The analysis of tree diversity revealed that disturbed forest plots have essentially the same tree species richness as undisturbed plots, 205 versus 206. Mean species richness per plot was 22 (SD 4) for the latter, 25 (SD 5) for the former, and not significantly different (t = −1.852, p value = 0.075). For large trees, however, undisturbed plots had 123 species, while disturbed plots only 94. With a mean species richness of 14 (SD 4) for undisturbed plots, and 11 (SD 3) for disturbed plots, the difference is significant (t = −2.745, p value = 0.010). In the case of small trees, across all five sites, disturbed plots house 144 species versus 128 in undisturbed plots, but the mean subplot species richness is 14 (SD 3) and 13 (SD 4), respectively, and not significantly different (t = 0.607, p value = 0.54).
The analysis of beta diversity revealed a high level of spatial dissimilarity across the five sites in Darien. Mean species turnover (Simpson dissimilarity) across all categories and permutations was 0.90 (min 0.78, max 0.98; SD 0.04), while nestedness was 0.01 (min 0.001; max 0.042; SD 0.006), thus species turnover is the primary driver of variation. Total beta diversity (Sorensen) was also very high with a mean of 0.91 (min 0.82; max 0.98; SD 0.032).

| The relative importance of Darien's AGB and tree diversity
We compared our plot-level AGB values with other wellstudied plots across the tropics ( Figure 5)  We also compared tree species richness in our plots in Darien with other plots across the Neotropics, (Table 2) (Condit, 1998;Hubbell et al., 1999Hubbell et al., , 2010Thompson et al., 2004;Valencia et al., 2004;Vallejo, Samper, Mendoza, & Otero, 2004). Noteworthy is the difference with the 50-ha BCI plot, which has 88 tree species ≥10 cm DBH fewer than our plots in Darien, despite having 22 additional ha (79% area increment). Our bootstrapped analysis revealed that BCI has 73-85 tree species (95% confidence intervals) ≥50 cm DBH per 30 randomly selected 1-ha plots, compared to169 species in our plots. In the case of subplots, BCI has 112-130 tree species (95% confidence intervals) ≥10 to <50 cm DBH, versus 219 in our subplots. It is important to note, however, that BCI is a single 50-ha plot located only a 1,560-ha island, while our 30 plots were distributed across a large landscape, thus likely capturing a greater proportion of species richness.

| Capturing AGB variation at the landscape level
Although the remote forests that we sampled had not undergone a significant and visible change in land use, such as to agriculture or cattle ranching, closer scrutiny revealed that not all these mature forests were intact. As mentioned above, half the plots had been subject to traditional indigenous extractive activities. We found that the key determinant of AGB variation is the level of disturbance. Neither differences across forest types, nor spatial variation, either random or due to biophysical factors, had as much effect on AGB as the selective extraction of large trees. Our results also highlight that this variation in ABG did not affect the forest canopy height, and that it was not captured by our satellitebased analyses of vegetation cover, or an analysis using LiDAR data (Asner et al., 2013).
Previous studies have highlighted that differences between disturbed and intact mature forests are not easily detected through remote sensing techniques unless the type, duration, frequency, intensity, and extent of disturbance results in significant changes to the forest canopy structure (Bustamante et al., 2016). Forest degradation indeed, unlike deforestation, occurs along a continuum (Putz & Redford, 2010), and its detection has proven challenging in the establishment of baselines for REDD+. Despite this, remote detection of degradation has been achieved in previous studies relying on time-series data from sites where changes in forest structure and cover are discernible, such as heavily logged and burned forests (Souza, Firestone, Silva, & Roberts, 2003), areas experiencing shifting cultivation (Pelletier, Codjia, & Potvin, 2012), or  where the presence of surface debris, bare soils and other human infrastructure (e.g., logging roads, skid tracks, adjacent agriculture) is detectable or used to infer degradation (Margono et al., 2012). This is not the case in the landscape that we studied where selective logging by Indigenous peoples does not lead to discernible damage to the forest or its surroundings. Therefore, disturbed forests may maintain the structure and characteristics of mature forests, even if AGB volumes differ vastly. The use of satellite time-series data, versus singletime frames has been proposed to identify more subtle changes (Pelletier et al., 2012), but our analyses of canopy height with two different time periods (i.e., 2005 and 2010) did not reveal significant differences. More extended time series might capture gap dynamics, or we may need to rely on proxies. For example, Bucki et al. (2012), in the context of forest monitoring for REDD+, proposed stratifying forested lands between intact and nonintact using distance from the forest edge as a proxy. In this case, the authors suggested a distance of 500 m. The distance proxy was based on results from previous studies that examined the penetration of impacts from various forms of forest disturbance.
In support of the use of proxies, we tested for and found significant differences between undisturbed and disturbed plots in terms of their distance to a nearest community (t = −2.0757, p value = 0.047). Undisturbed plots were on average at 4.6 km from the nearest community, while at 3.1 km for disturbed forest plots. This is consistent with previous studies in eastern Panama that examined the relationship between the harvest of useful plants and distance to a community (Dalle & Potvin, 2004). Thus, for remote but inhabited mature forests, one could conceivably define a "buffer zone" of disturbed forests around settlements, and consider that the remaining forest matrix has higher carbon stocks than areas surrounding communities.

| Climate change mitigation and biodiversity conservation in Darien
At the global scale, carbon-rich and high-biodiversity areas overlap (Strassburg et al., 2010;Venter et al., 2009), but the trend is not universal (Murray, Grenyer, Wunder, Raes, & Jones, 2015). Our study revealed that even when disturbed forest plots have up to 54% less carbon than undisturbed plots, they maintain the same tree species richness. As such, policy instruments with a focus on carbon or biodiversity, may not be able to fully capitalize on "filling two needs with one deed." Instead, policy development will have to explicitly factor both biodiversity conservation and carbon sequestration, since an umbrella effect in which one can be a surrogate for the other may not always exist. It is important to note, however, that even though disturbed forest areas may not hold the same volume of forest carbon stocks, studies (Bongers, Chazdon, Poorter, & Peña-Claros, 2015;Chazdon et al., 2016) have revealed their disproportionately high capacity to sequester CO 2 . Thus, a lack of overlap between carbon and species rich areas does not necessarily preclude these disturbed forest areas from being targeted by policy instruments such as REDD+.
Darien was identified as a priority area during Panama's REDD-Readiness process (ANAM, 2008). But how REDD+ may be operationalized in this region remains an open question. Our results position Darien's forests among the most carbon-rich in the tropics, and most tree species rich in the Neotropics. These findings underscore the global importance of this region for climate change mitigation and conservation. However, Darien is undergoing a process of rapid transformation, colonization and deforestation (Arcia Jaramillo, 2015). The region is considered among the most threatened frontier forests in the Americas (Bryant et al., 1997;WWF, 2015) due to century-old plans to open the "Darien Gap" (Miller, 2014), and more recently, the Panamanian and Colombian governments' plans to develop an electrical interconnection project that could pass through Darien (Proyecto Mesoamérica, 2015). Any linear infrastructure, whether roads or power lines, could lead to further colonization and exploitation (Laurance, Goosem, & Laurance, 2009).
In regions inhabited by Indigenous peoples or local groups, clear and secure land tenure has been posited as a key requirement for REDD+ (Larson et al., 2013). In the case of Darien, most forested areas are located on lands under Indigenous stewardship (Vergara-Asenjo & Potvin, 2014), many that overlap fully or partially with protected areas under different management categories. This complex matrix of overlapping land tenure and management regimes, coupled with a long history of land invasions by migrant colonist farmers (Herlihy, 2003), creates a challenging sociopolitical context for REDD+, but one that the proponents of the mechanism must address if they intend to deliver on expected REDD+ outcomes and cobenefits.
The importance and value of intact and near intact forests, those that have been mostly spared from anthropogenic impacts, has been underscored in previous studies (Potapov et al., 2017). These forests, including Darien, make significant contributions, in comparison to degraded forests, to conserving biodiversity, sequestering and storing carbon, and providing a variety of other key ecosystem services. In addition to these global environmental values, intact forest areas may also hold fundamental cultural and social values, particularly for Indigenous peoples (Watson et al., 2018). Darien's three indigenous cultures, namely the Emberá, Wounaan and Guna, are intimately tied to these forests (Heckadon-Moreno, Herrera, & Pastor Nuñez, 1984). This study provides supporting evidence for the imperative need to direct greater attention to the intact forests of Darien, given the disproportionately important role they play in mitigating climate change and housing biodiversity. Darien's undisturbed forests have the highest carbon stocks among nine mature forest sites across the Neotropics, and the second highest tree species richness among five mature forest sites in the same region.