**Assessing Tree Coverage and the Direct and Mediation Effect of Tree Diversity on Carbon Storage through Stand Structure in Homegardens of Southwestern Bangladesh**

**Md Mizanur Rahman 1,2,† , Gauranga Kumar Kundu 3,†, Md Enamul Kabir <sup>3</sup> , Heera Ahmed <sup>3</sup> and Ming Xu 1,2,\***


**Abstract:** Dealing with two major challenges, climate change mitigation and biodiversity loss, under the same management program, is more noteworthy than addressing these two separately. Homegardens, a sustainable agroforestry system and a home of diverse species, can be a possible choice to address these two issues. In this study, we assessed tree coverage, and the direct and indirect effects of tree diversity on carbon storage in different carbon pools through stand structure in homegardens of southwestern Bangladesh, using Sentinel 2 and field inventory data from 40 homesteads in eight villages. An unsupervised classification method was followed to assess homegardens' tree coverage. We found a high tree coverage (24.34% of total area of Dighalia) in homesteads, with a high overall accuracy of 96.52%. The biomass and soil organic carbon (*p* < 0.05) varied significantly among the eight villages, while total carbon stock did not vary significantly (*p* > 0.05). Shannon diversity had both direct and indirect effects on biomass carbon, upper layer soil organic carbon and total carbon storage, while basal area mediated the indirect effect. Both basal area and tree height had positive effects on biomass carbon and total carbon storage, with basal area having the strongest effect. These findings suggest that we must maintain higher diversity and tree height in order to maximize and sustain carbon storage, where tree diversity increases stand basal area and improves total carbon storage (including soil organic) in homegardens. Therefore, privately managed homegardens could be a potential nature-based solution for biodiversity conservation and climate change mitigation in Bangladesh.

**Keywords:** traditional homegardens; agroforestry system; biodiversity and carbon; optical remote sensing

#### **1. Introduction**

Nature-based solutions are one of the best strategies for adapting to global climate change and biodiversity loss [1]. To combat global warming, world leaders committed in the Paris Agreement to keep atmospheric temperature below 2 ◦C by end of the century [2]. Each signatory country to the Paris Agreement has pledged to cut a part of its greenhouse gas emissions by 2030 [2]. Under the framework of Nationally Determined Contributions (NDCs), each country is meeting their obligation by lowering either low carbon sectors (transport, industry, power, etc.) or by increasing forestry activities (afforestation or restoration) [3]. Homegardens, a sustainable and integrated agroforestry system in tropical and subtropical countries, has huge carbon sequestration potential [4,5]. This climate regulatory role of homegardens can be a potential NDC component for many countries in meeting its carbon emission commitment [3].

**Citation:** Rahman, M.M.; Kundu, G.K.; Kabir, M.E.; Ahmed, H.; Xu, M. Assessing Tree Coverage and the Direct and Mediation Effect of Tree Diversity on Carbon Storage through Stand Structure in Homegardens of Southwestern Bangladesh. *Forests* **2021**, *12*, 1661. https://doi.org/ 10.3390/f12121661

Academic Editors: Himlal Baral and Mi Sun Park

Received: 1 October 2021 Accepted: 23 November 2021 Published: 30 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Sustainably managed lands that have the ability to supply food and other domestic products to conserve biodiversity and regulate climate are of high demand in the current changing climate [6]. Homegardens are a tiny piece of land just in and around the homestead in rural areas of the tropical and subtropical regions. This unique type of agroforestry system provides multiple ecosystem services over the year to members of the household [7]. It is an intensive cropping system of domesticated annual and perennial plants and/or animals, which are the primary source of household demand, such as food, fruits, fodder, fuel wood and furniture [7,8]. Different coexisting species of trees, shrubs and herbs in homegardens with vertical differentiation make this traditional agroforestry system a forest-like structure [7]. Within human-dominated landscapes, homegardens are considered as highly biodiverse agroforestry systems which play an important role in biodiversity conservation and carbon sequestration [4].

Studies on homegardens' carbon storage have shown a high capacity for minimizing atmospheric carbon through photosynthesis compared to other agroforestry systems [4,9,10]. Findings on homegardens carbon stocks have revealed that aboveground carbon is mainly contributed by tree and palm species in homegardens, and that soil carbon occupies the largest percentage of the total carbon stored in homegardens [9]. Homegardens with high species diversity and stem density lead to a more basal area and biomass (above and below ground), and therefore, they store more biomass and soil carbon compared to homegardens with less species density and stem density [5,11,12]. In the forest ecosystem, the positive role of species diversity in carbon storage is underpinned by the niche complementarity hypothesis, which believes that higher numbers of species in a stand with greater trait variation have used more resources, and therefore, store more carbon [13]. Some studies have, however, documented no particular relationship between species diversity and carbon stock, which can be explained by selection effect. The selection effect believes that among the coexisting species, the most dominant species with its key traits such as size vigor (height and diameter) regulates carbon storage [12]. Both of these hypotheses can be a driver for homegarden agroforestry carbon regulation, because it is a multi-storiedspecies landscape. Furthermore, inconsistent results of species diversity and carbon storage relationships have also been reported in different studies in agroforestry systems [5,14]. Thus, for understanding the species diversity and carbon storage and their relationship in a specific study area, we need to conduct direct field investigations.

Bangladesh committed to the Paris Agreement to reduce its greenhouse gas emission from power, transport and industry sectors with an equivalent to 12 Mt CO2e by 2030, which is 5% below the total business as usual emission level from those sectors. However, Bangladesh is a country of villages where almost every household has a homegarden, which are well-established land use systems [7]. A recent study reported that the total tree coverage in non-forest area is about 4.5 million ha, where homegardens have dominant contributions [15]. The homegardens in Bangladesh are biologically diverse and most of the studies on homegardens in Bangladesh remained descriptive on the floristic characteristics, structure, uses and the relation between household and homegarden characters [7,16]. Homegardens in Bangladesh can be a potential component of NDCs, similar to other developing countries. However, few studies focused on the total carbon storage capacity and carbon and biodiversity relationship in Bangladesh [17,18]. This knowledge gap in homegardens agroforestry system in Bangladesh might be one of the most important reasons for excluding agroforestry as a key NDCs component by policy makers, despite its huge coverage and high potential for climate change mitigation.

In this study, we assessed homegarden tree coverage, stand structure, species diversity and carbon storage in southwestern Bangladesh. We set three specific objectives in this study. First, to assess the total tree coverage in homesteads of Dighalia upazila (administrative unit) in the Khulna district, Bangladesh by using Sentinel 2 imagery. Second, to quantify the direct and indirect effects of species diversity on biomass and soil organic carbon stocks through different stand structures in homegardens. Third, to ex-

amine how results in line with carbon emission targets have contributed to homegardens' management system. results in line with carbon emission targets have contributed to homegardens' management system.

to quantify the direct and indirect effects of species diversity on biomass and soil organic carbon stocks through different stand structures in homegardens. Third, to examine how

*Forests* **2021**, *12*, x FOR PEER REVIEW 3 of 16

#### **2. Materials and Methods 2. Materials and Methods**

#### *2.1. Study Site 2.1. Study Site*

The study was conducted in Dighalia Upazilla (administrative unit) in the Khulna district of Bangladesh (Figure 1). Dighalia Upazilla is located in the southwestern part of Bangladesh. The study area is primarily a flood plain land mass between 22.50◦ to 22.60◦ N latitude and 89.31◦ to 89.37◦ E longitude. The deltaic landscape of this region is mainly low (<10 m above average sea level), flat and fertile [19]. The average area of the homegarden over the study site is 0.05 ha and rich in biodiversity [19]. It generally enjoys a tropical to subtropical monsoon climate with an average annual temperature of 26 ◦C. January is the coolest month (average temperature: 12.4 ◦C), while April is the hottest month (average temperature: 34.6 ◦C) in this region. The average annual rainfall for this region is 1986 mm (range: 1400–2600 mm). This upazila is enclosed by three main rivers, i.e., the Bhairab, Citra and Naboganga, and is adjacent to Khulna City, the largest city in southwestern Bangladesh [20]. The study was conducted in Dighalia Upazilla (administrative unit) in the Khulna district of Bangladesh (Figure 1). Dighalia Upazilla is located in the southwestern part of Bangladesh. The study area is primarily a flood plain land mass between 22.50° to 22.60° N latitude and 89.31° to 89.37° E longitude. The deltaic landscape of this region is mainly low (<10 m above average sea level), flat and fertile [19]. The average area of the homegarden over the study site is 0.05 ha and rich in biodiversity [19]. It generally enjoys a tropical to subtropical monsoon climate with an average annual temperature of 26 °C. January is the coolest month (average temperature: 12.4 °C), while April is the hottest month (average temperature: 34.6 °C) in this region. The average annual rainfall for this region is 1986 mm (range: 1400–2600 mm). This upazila is enclosed by three main rivers, i.e., the Bhairab, Citra and Naboganga, and is adjacent to Khulna City, the largest city in southwestern Bangladesh [20].

**Figure 1.** Map showing carbon sampling plots and ground truthing points in Dighalia Upazilla.

#### *2.2. Sampling Design and Field Inventory*

A multi-stage random sampling procedure was followed to select homegardens. First, two unions (smallest administrative unit), namely Senhati and Barakpur, from the total of six unions of Dighalia Upazila were selected randomly. Then, eight villages were selected randomly from the two randomly selected unions (four villages from each union). A total of 40 plots of 10 × 10 m (five from each randomly selected village) were taken purposively (homestead with homegarden), as shown in Figure 1. We measured the diameter at breast height (DBH) and canopy height of all species with a DBH ≥ 3 cm. All species in each sampled plot were identified and recorded to species level or by local name and were later confirmed from an authentic source [19]. We collected soil samples for soil organic carbon assessment from the center of each plot. We pulled out soil cores using an open face peat augur. From the soil core, we took a total of 80 (40 homegardens × 1–plot × 2 depths) soil samples; two from each plot using a five-centimeter-long steel core at the midpoint of 0–15 cm depth and 15–30 cm depth. All samples were air-dried and analyzed for estimating bulk density and soil organic carbon. We analyzed all the soil samples at the Nutrient Dynamics Laboratory of Khulna University, Bangladesh.

#### *2.3. Remote Sensing Data and Processing*

Sentinel -2 satellite image that was acquired in 13 December 2016 was used in this study for assessing the vegetation coverage. The Level-1C product of Sentinel 2 is a version of Top-Of-Atmosphere reflectance image, which is radiometrically and geometrically corrected. Sentinel 2 is a multispectral imager launched on 23 June 2015 by the European Space Agency, which has 13 bands of different resolutions (10–60 m) [21]. In this study, we used the 10 m resolution bands (Red and Near Infrared) for the Normalized Difference Vegetation Index (NDVI). We masked non-homestead areas (water body, crop field and fishpond) in order to distinguish rural settlement with homesteads. Then, we applied an ISODATA algorithm, an iteration method that collapses clusters of similar groups into one by measuring the Euclidean distance between the cluster centers [22,23]. We set the number of classes in the classification to 10, the total class size to 10 and the sample interval to 2. Using Google Earth visual interpretation and field data, we assigned homegardens and non-homegardens classes (houses, paddy fields, fish ponds, waterbodies and built-up areas) to the resultant one band image of ISODATA classification with 10 clusters [24]. We collected 115 points randomly from homegardens (59) and non-homegardens (56) within homestead areas (Figure 1) by using Google Earth visual interpretation for accuracy assessment of our classification result [23–25]. A confusion matrix was used for assessing the accuracy of our classification [25]. We used ENVI 5.00 and QGIS 3.14 for classification and mapping.

#### *2.4. Data Analysis*

#### 2.4.1. Biomass Carbon

Common allometric equations were used (Table 1) to calculate the biomass [26,27], as it was developed from wide graphical (tropical countries) data and has diameter range (Table 1). The wood density data were obtained from the Global Wood Density Database [28]. For species with missing wood density in this database, we used the average wood density. Biomass carbon was calculated by multiplying 0.5, as it is assumed that wood biomass contains 50% carbon [28].

**Table 1.** Allometric equations used for biomass estimation in this study.


AGB = Aboveground biomass, ρ = Wood density (gcm−<sup>3</sup> ), DBH = Diameter at breast height, ht = Height, BGB = Belowground Biomass.

#### 2.4.2. Soil Carbon Calculation

Soil carbon storage was estimated by the product of soil organic carbon concentration, soil bulk density and soil depth range. Soil bulk density was determined for each soil layer by dividing the oven-dried soil sample mass by the volume of the sample (Equation (1)) [30]:

 $Bulk\ Density = \frac{\left(Wt\_{105^{\circ}C}\right)}{V\_{\text{core}}}$ 
$$V\_{\text{core}} = \frac{\left(\pi D\_{\text{core}}2L\_{\text{core}}\right)}{4}$$

where *Wt*105◦<sup>C</sup> is the weight of oven dried soil, *Vcore* is the volume of the core, *Dcore* is the inner diameter of the core and *Lcore* is the length of the core.

The Loss of Ignition method was used for calculating soil organic matter [30]. Generally, soil organic matter contains 58% organic carbon [31]. Thus, we divided the organic matter by 1.724 (a universal conversion factor (called the van Bemmelen factor) for converting the soil organic carbon percentage [31].

After getting the soil bulk density and organic carbon percentage, we calculated the soil carbon stock for each layer using Equation (2):

$$\text{Soil C (Mg ha}^{-1}\text{)}=\text{bulk density (g m}^{-3}\text{)} \times \text{soil depth interval (m)} \times \% \text{OC} \times 0.01 \tag{2}$$

where Soil depth interval = 0.15 for 0–15 cm depth and 15–30 cm depth, %OC is expressed as a decimal fraction (e.g., 5% is expressed as 0.05) and 0.01 is a conversion factor to convert units to Mg ha−<sup>1</sup> .

#### 2.4.3. Total Carbon Stock

We calculated total carbon stocks (Equation (3)) per hectare by summing the biomass carbon stocks (aboveground and belowground biomass carbon) and soil organic carbon stocks (soil carbon, 0–15 cm and soil carbon, 15–30 cm):

Total C stock (Mg ha−<sup>1</sup> ) = Biomass carbon(AGC + BGC) + TSOCSOC (0–15cm) + SOC (15–30) (3)

#### 2.4.4. Woody Species Diversity Calculation

The Shannon–Wiener index (*H*; Equation (4)) was used for species diversity [32], while for species richness, the Margalef index (*D*; Equation (5)) was used [33]:

$$H = \sum\_{i=1}^{S} (P\_i \times \ln(P\_i))\tag{4}$$

$$D = \frac{(\mathbb{S} - 1)}{\ln(N)} \tag{5}$$

where *H* is the Shannon diversity for a plot, *S* is the number of species (species richness), *Pi* is the proportional of individuals of species *i* in the plot, *N* is the total number of individuals in a plot and "ln" is the natural logarithm.

One-way analysis of variance (ANOVA) was used to test the significant difference of carbon pools across the eight villages and Duncan's Multiple Range Test was used for multiple comparisons. Pearson's correlation tests were performed using the "rstatix" package in R environment (Version 3.6.3) [34], to explore the correlation between stand structure, woody species diversity and carbon pools. Structural equation model (SEM) was used to quantify the direct and indirect effects of tree diversity on carbon pools through stand structure [35]. We used the Lavaan package for employing SEM in R environment (Version 3.6.3) [34]. As recommended, we log transformed and standardized all the covariables before applying the structural equation model. Multiple linear regression was used to quantify the variance inflation factors (VIF) with a threshold value of 2 and assuming that covariables having a value below this threshold have no multicollinearity [36]. In the final

steps of the variance inflation factors test, mean tree height, basal area, stem density and Shannon diversity were retained (MDBH had a VIF value that exceeded the threshold value, Supplementary Table S1). Then, we employed SEM for quantifying the direct, indirect and total effect of the Shannon diversity on different carbon pools through stem density, canopy height and basal area. A total 14 SEMs were tested for seven carbon pools, where highly insignificant covariables were removed until the model satisfied fit statistics [37] (Supplementary Table S2). The indirect effects of Shannon diversity were assessed by the path value through its mediator to different carbon pools (response variable) [35]. We assessed the total effect of the Shannon diversity on a carbon pool by adding the direct and indirect effect through different mediators [35,37].

#### **3. Results**

#### *3.1. Homegardens Cover*

The spatial distribution of homegardens and non-homegardens based on the unsupervised classification of Sentinel 2 images in Dighalia Upazilla showed a good match with the Google Earth map (same year image as of Sentinel 2; Figure 2). Our classification separated homegardens and non-homegardens in homestead areas in Dighalia with a higher overall accuracy and kappa coefficient (overall accuracy = 96.52% and Kappa = 0.93; Table 2). The overall coverage of the homegardens was 17.32 km<sup>2</sup> and the total percentage of homegardens in Dighalia was 24.34 (the total area of Dighalia is 71.16 km<sup>2</sup> ). *Forests* **2021**, *12*, x FOR PEER REVIEW 7 of 16

**Figure 2.** Extent of homegardens and non-homegardens in Dighalia Upazila. (**a**) Based on unsupervised classification of Sentinel 2 imagery; (**b**) Google Earth satellite image® shows homegardens covered by dark area. **Figure 2.** Extent of homegardens and non-homegardens in Dighalia Upazila. (**a**) Based on unsupervised classification of Sentinel 2 imagery; (**b**) Google Earth satellite image® shows homegardens covered by dark area.

**Table 2.** Result of the classification accuracy based on the Confusion Matrix. Column represents the Google Earth vegetation class, while the row represents the Sentinel 2-based class. **Table 2.** Result of the classification accuracy based on the Confusion Matrix. Column represents the Google Earth vegetation class, while the row represents the Sentinel 2-based class.


13.76 ± 0.96 m2 ha−1, respectively, across the 40 sample plots (Table 3).

*3.2. Diversity and Stand Structure* 

Producers (%) 94.92 98.21 Overall 96.52

The mean Shannon diversity, Margalef diversity and stand structures exhibited a spatial variation across the 40 plots (Table 3). The range of Shannon diversity was 1.91 to 3.01 with a mean of 2.56 ± 0.04 in the 40 plots (Table 3). The Margalef diversity index had a mean value of 2.31 ± 06 (range: 1.44–3.03). The estimated mean stem density, DBH, tree height and basal area were 1055.00 ± 35.62 treesha−1, 13.30 ± 0.49 cm, 7.88 ± 0.17 m and

#### *3.2. Diversity and Stand Structure*

The mean Shannon diversity, Margalef diversity and stand structures exhibited a spatial variation across the 40 plots (Table 3). The range of Shannon diversity was 1.91 to 3.01 with a mean of 2.56 ± 0.04 in the 40 plots (Table 3). The Margalef diversity index had a mean value of 2.31 ± 06 (range: 1.44–3.03). The estimated mean stem density, DBH, tree height and basal area were 1055.00 <sup>±</sup> 35.62 treesha−<sup>1</sup> , 13.30 ± 0.49 cm, 7.88 ± 0.17 m and 13.76 <sup>±</sup> 0.96 m<sup>2</sup> ha−<sup>1</sup> , respectively, across the 40 sample plots (Table 3).

**Table 3.** Descriptive statistics of stand structures and diversity indices across the 40 plots in Dighalia.


#### *3.3. Carbon Stocks*

High variations of above- and belowground biomass carbon were observed (aboveground: 45.01 <sup>±</sup> 6.95 to 88.00 <sup>±</sup> 6.39 Mg ha−<sup>1</sup> ; belowground: 9.98 <sup>±</sup> 1.36 to 18.11 <sup>±</sup> 2.6 Mg ha−<sup>1</sup> ). We found similar biomass carbon stocks in Senhati, Baracpur and Lakhoati villages, but these three villages have significantly lower biomass carbon than in the Chandani mahal village (*p* < 0.05; Table 4). However, the biomass carbon stocks in these four villages were not significantly different from the other four villages (*p* > 0.05). Overall, aboveand belowground carbon stocks across the 40 plots in Dighalia were 60.42 ± 3.93 and 12.89 <sup>±</sup> 0.75 Mg ha−<sup>1</sup> , respectively. Although the upper layer (0–15 cm) soil organic carbon stocks were insignificant across the different villages (*p* > 0.05; Table 4), the lower layer (15–30 cm) soil organic carbon stocks were found to be significantly different (*p* < 0.05; Table 4). In the Chandani mahal and Hagigram villages, soil organic carbon stocks at lower layers were significantly lower than at the Bativita village (*p* < 0.05; Table 4). Although we found a significant difference in biomass and lower layer soil organic carbon stocks, the differences in overall carbon stocks were insignificant across the villages (*p* > 0.05). The total carbon stock (biomass and soil carbon pools) in the eight villages ranged from 109.10 <sup>±</sup> 10.89 Mg ha−<sup>1</sup> (in Senhati) to 152.22 <sup>±</sup> 10.15 Mg ha−<sup>1</sup> (in Chandani mahal) with an overall average total carbon stock of 129.47 <sup>±</sup> 5.10 Mg ha−<sup>1</sup> (Table 4). The overall biomass carbon stock (73.31 <sup>±</sup> 4.67 Mg ha−<sup>1</sup> ) contributes 56.62 percent of the total carbon stock, while soil carbon (56.16 <sup>±</sup> 1.71 Mg ha−<sup>1</sup> ) contributes 43.38 percent of the total carbon stock. As we mapped the total area of homegardens and mean total carbon per hectare in Dighalia, we also assessed the total carbon stock in Dighalia by multiplying the mean carbon stocks by the total area. The total carbon stocks in Dighalia homegardens was 224,242.2 Mg C, which is equivalent to 822,968.29 Mg CO2.

**Table 4.** Mean (±SE) difference of different carbon pools between villages. According to the Duncan's Multiple Range Test, the columns with a similar letter suggest a non-significant difference (*p* > 0.05).


#### *3.4. Bivariate Relationship between Stand Structure, Diversity and Carbon Pools* (aboveground and belowground), where basal area had a stronger relation compared to

*3.4. Bivariate Relationship between Stand Structure, Diversity and Carbon Pools* 

*Forests* **2021**, *12*, x FOR PEER REVIEW 9 of 16

All four stand structure parameters were positively related to biomass carbon (aboveground and belowground), where basal area had a stronger relation compared to stem density, mean tree height and mean diameter (*p* < 0.05; Figure 3). However, these four-stand structure parameters had no significant relation with any of the soil organic carbon layers (*p* > 0.05; Figure 3). When we considered all the carbon pools together as total mean carbon stocks, all four stand structure parameters except mean tree height had a significant relation and basal area also had a stronger relation (*p* < 0.05, Figure 3). Shannon diversity index also had positive effects on all carbon pools (except lower layer soil organic carbon; *p* > 0.05; Figure 3) and total carbon stocks as well as stem density and basal area (*p* < 0.05; Figure 3). stem density, mean tree height and mean diameter (*p* < 0.05; Figure 3). However, these four-stand structure parameters had no significant relation with any of the soil organic carbon layers (*p* > 0.05; Figure 3). When we considered all the carbon pools together as total mean carbon stocks, all four stand structure parameters except mean tree height had a significant relation and basal area also had a stronger relation (*p* < 0.05, Figure 3). Shannon diversity index also had positive effects on all carbon pools (except lower layer soil organic carbon; *p* > 0.05; Figure 3) and total carbon stocks as well as stem density and basal area (*p* < 0.05; Figure 3).

All four stand structure parameters were positively related to biomass carbon

**Figure 3.** Correlation matrix between stand structure, diversity and different carbon pools. Block with cross sign indicates insignificant relation (*p* > 0.05). MDBH: Mean Diameter at Breast Height (1.37 m), MTH: Mean tree height, SD: Stem Density, BA: Basal Area, D = Margalef Diversity index, H = Shannon Diversity Index, AGC: Aboveground Carbon, BGC: Belowground Carbon, SOC1 = Soil Organic Carbon at 0–15 cm depth, SOC2 = Soil Organic Carbon at 15–30 cm, TBC = Total **Figure 3.** Correlation matrix between stand structure, diversity and different carbon pools. Block with cross sign indicates insignificant relation (*p* > 0.05). MDBH: Mean Diameter at Breast Height (1.37 m), MTH: Mean tree height, SD: Stem Density, BA: Basal Area, D = Margalef Diversity index, H = Shannon Diversity Index, AGC: Aboveground Carbon, BGC: Belowground Carbon, SOC1 = Soil Organic Carbon at 0–15 cm depth, SOC2 = Soil Organic Carbon at 15–30 cm, TBC = Total Biomass Carbon (AGC + BGC), TSOC = Total Soil Organic Carbon (0–30 cm).

Biomass Carbon (AGC + BGC), TSOC = Total Soil Organic Carbon (0–30 cm).

#### *3.5. Direct, Indirect and Total Effects of Biodiversity and Stand Structure on Carbon Pools 3.5. Direct, Indirect and Total Effects of Biodiversity and Stand Structure on Carbon Pools*

Shannon diversity had significant direct and indirect effects on aboveground biomass carbon, belowground biomass carbon total biomass carbon, and total carbon storage, but had only direct effects on the upper layer soil organic carbon and total soil organic carbon (Figure 4a–g). The basal area only significantly mediated the Shannon diversity effect on biomass carbon (above- and belowground) and total carbon stocks (Figure 4a,b,e,g). Basal area had the strongest direct effect on biomass and total carbon storage (Figure 4a,b,e,g), while Shannon diversity had the strongest effect on upper layer and total soil organic carbon storage (Figure 4c,f). Tree height also had positive effects on biomass carbon and total carbon storage (Figure 4a,b,e,g). As the Shannon diversity has indirect effects on the biomass and total carbon storage and, which were mediated by the basal area, the total effect of the Shannon diversity on the aboveground biomass carbon, belowground biomass carbon, total biomass carbon and total carbon storage were 0.393, 0.393, 0.513, 0.393, 0.411 and 0.544, respectively (Table 5). The Shannon diversity, basal area and mean tree height together explained 80%, 80%, 27%, 04%, 80%, 19% and 61% ariations in aboveground biomass carbon, belowground biomass carbon, upper layer soil organic carbon, total biomass carbon, total soil organic carbon and total carbon storage in homegardens, respectively (Figure 4a–g). Shannon diversity had significant direct and indirect effects on aboveground biomass carbon, belowground biomass carbon total biomass carbon, and total carbon storage, but had only direct effects on the upper layer soil organic carbon and total soil organic carbon (Figure 4a–g). The basal area only significantly mediated the Shannon diversity effect on biomass carbon (above- and belowground) and total carbon stocks (Figure 4a,b,e,g). Basal area had the strongest direct effect on biomass and total carbon storage (Figure 4a,b,e,g), while Shannon diversity had the strongest effect on upper layer and total soil organic carbon storage (Figure 4c,f). Tree height also had positive effects on biomass carbon and total carbon storage (Figure 4a,b,e,g). As the Shannon diversity has indirect effects on the biomass and total carbon storage and, which were mediated by the basal area, the total effect of the Shannon diversity on the aboveground biomass carbon, belowground biomass carbon, total biomass carbon and total carbon storage were 0.393, 0.393, 0.513, 0.393, 0.411 and 0.544, respectively (Table 5). The Shannon diversity, basal area and mean tree height together explained 80%, 80%, 27%, 04%, 80%, 19% and 61% ariations in aboveground biomass carbon, belowground biomass carbon, upper layer soil organic carbon, total biomass carbon, total soil organic carbon and total carbon storage in homegardens, respectively (Figure 4a–g).

**Figure 4.** Structural equation models. (**a**–**g**) Aboveground biomass carbon (**a**), belowground biomass carbon (**b**), soil organic carbon (0–15 cm) (**c**), soil organic carbon (15–30 cm) (**d**), total biomass carbon (**e**), total soil organic carbon (**f**) and total carbon (**g**). All SEMs had a similar insignificant χ2 (Chi-Square) of 1.22 (*p* = 0.268), with a comparative fit index close to one (CFI; 0.98–0.99) and a standardized root mean square residual close to zero (0.08), indicating no significant deviation from model and datasets at one degree of freedom. The lines with pink and black indicate a negative and positive association between the two covariables. Arrows with numbers indicate the standardized association of predictors with dependent variables. Numbers with percentages above boxes independent variables indicate their explained variance (Coefficient of determinant: R squared) by all the predictors. The path values with asterisks indicate their significance level (\*\*\* *p* < 0.001, \*\* *p* < 0.01, \* *p* < 0.05), while the insignificant paths were indicated as dotted lines (*p* > 0.05). **Figure 4.** Structural equation models. (**a**–**g**) Aboveground biomass carbon (**a**), belowground biomass carbon (**b**), soil organic carbon (0–15 cm) (**c**), soil organic carbon (15–30 cm) (**d**), total biomass carbon (**e**), total soil organic carbon (**f**) and total carbon (**g**). All SEMs had a similar insignificant χ 2 (Chi-Square) of 1.22 (*p* = 0.268), with a comparative fit index close to one (CFI; 0.98–0.99) and a standardized root mean square residual close to zero (0.08), indicating no significant deviation from model and datasets at one degree of freedom. The lines with pink and black indicate a negative and positive association between the two covariables. Arrows with numbers indicate the standardized association of predictors with dependent variables. Numbers with percentages above boxes independent variables indicate their explained variance (Coefficient of determinant: R squared) by all the predictors. The path values with asterisks indicate their significance level (\*\*\* *p* < 0.001, \*\* *p* < 0.01, \* *p* < 0.05), while the insignificant paths were indicated as dotted lines (*p* > 0.05).

**Table 5.** Indirect and total standardized effect of the Shannon diversity on different carbon pools. The path values with asterisks indicate their significance level (\*\*\* *p* < 0.001, \*\* *p* < 0.01, \* *p* < 0.05).


#### **4. Discussion**

Homegardens' agroforestry system is rich in biodiversity in the tropics and provides a variety of ecosystem services, including carbon sequestration. We quantified homegardens' agroforestry vegetation coverage and the direct and mediation effects of tree diversity on carbon storage via stand structure in Dighalia, Khulna, Bangladesh. We found that homegardens had high vegetation coverage and carbon stocks where tree diversity promotes carbon storage directly and indirectly through basal area. The stand basal area had a stronger effect on carbon storage than other stand structural parameters.

#### *4.1. Homegardens Vegetation Coverage*

We found a high vegetation coverage (as homegardens) in homesteads, which accounted for 24% of the total area of Dighalia Upazila. Two recent studies reported that the percentage of tree coverage outside forests (woodlots, homegardens and other plantations) in Bangladesh is 13.01 (1,920,700 ha/country area) by Landsat Satellite with 30 m × 30 m pixel [38] and 15.1% by Sentinel 1 and 2 with 10 m × 10 m pixel [15]). While the other two studies focus on all types of trees outside forests, we focused only on the homegardens that were found in villages, which play a key role in biodiversity conservation and carbon sequestration. The difference in percentage of tree coverage in our study could be due to the spatial resolution of satellite data and the specific types of trees outside forests. The forest area percentage within the total land area is one of the key indicators of the Sustainable Development Goals (15.1.1), Achi biodiversity targets of the convention on biological diversity (Target 5; related to forest coverage or extent) and a key activity data for estimating countries' emission factor [6]. Higher accuracy in the estimation of homegardens' coverage with high-resolution imagery (Sentinel 2 10 m band), thus, contributes to assessing the accurate uncertainty in the carbon emission factor.

#### *4.2. Stand Structure*

Forest stand structures (tree height, DBH, stand density and basal area) are the key indicators of growth and productivity modeling, and recently, the stand level means of these parameters have been used in stand level biomass carbon modeling [39]. Higher amounts of those parameters at the stand level reveals standing resources and, thereby, growth and biomass carbon stocks [18,39]. The mean stem density in our study (1055.00 <sup>±</sup> 35.62 trees ha−<sup>1</sup> ) was very high compared with the natural forests in the Chittagong Hill Tracts (South) Forest Division (381 trees ha−<sup>1</sup> ; Nath et al., 1998) in the Chunati Wildlife Sanctuary, Cox's Bazar (459 trees ha−<sup>1</sup> ) [40]. Higher stem density is also reported in homegardens in Ethiopia (1125.23 ± 334.6) [41]. Homegardens are an intensive land use system, where low planting space is being maintained, which may be the cause of the higher stem density in homegardens compared to natural and secondary forests in Bangladesh. These traditional homegardens with high stem density, thus, play a major role in meeting the timber demand in a sustainable manner while reducing pressure on natural forests [7,18].

#### *4.3. Carbon Stocks*

We found a high carbon stock in the studied homegardens, which indicates a huge potential for homegardens to mitigate climate change through sequestering atmospheric CO2. The total average carbon stock in our study (129.47 <sup>±</sup> 5.10 Mg ha−<sup>1</sup> ; biomass + soil 0–30 cm) was within the range of total carbon stocks in Ethiopian homegardens (91.75 ± 4.31 to 156.17 <sup>±</sup> 13.78 Mg ha−<sup>1</sup> ) [14], but higher than Indonesian homegardens (107 Mg ha−<sup>1</sup> ) [9]. Our result was also higher than that of other agroforestry system's carbon stocks, such as agroforest (120 Mg ha−<sup>1</sup> ; soil carbon at 40 cm only) in Panama and traditional agroforestry systems (14–70.08 Mg ha−<sup>1</sup> ) in Mali [42]. The variations in total carbon stocks in different studies are mainly due to the difference in biomass and soil carbon pools that depend on higher stem density, basal area and tree size (height; Table 4 and Figure 3) as well as homestead age [9,11]. For example, in our study, the stem density (1055 stem/ha) in homegardens was higher than that of Indonesian homegardens (624.4 stem/ha). Our findings also showed that across villages, some had higher biomass carbon and lower soil carbon stocks than other villages and vice versa, where the biomass carbon varies mainly due to the difference in mean stem density, basal area and canopy height (Table 3; Figure 3). The positive correlation between these stand structures and different carbon pools indicates that high carbon storage is present in homegardens, and thus, villages are associated with higher stand structure (Figure 3).

#### *4.4. Effects of Tree Diversity and Stand Structure on Carbon Pools*

Dealing with two major challenges, climate change mitigation and biodiversity loss under the same management program is more noteworthy than addressing these two separately [43]. In order to practice this type of management, we need to show the evidence that species diversity and carbon storage are positively linked. However, conflicting results (positive, negative or neutral) in different ecosystems (or forest types), sparked debate among scientists and posed a barrier to the adoption of a unified management plan [5,11,18,44]. In our study, we found that tree diversity (Shannon) promotes all carbon pools (except for lower layer soil organic carbon) including total carbon stock (Figures 3 and 4). Tree diversity also had a mediation effect on all carbon pools (except soil organic carbon) through increasing stand basal area. These findings indicate that homegardens associated with high diversity had higher carbon stocks. Our result, therefore, supports the complimentary niche hypothesis which believes that higher niche differences between coexisting species can promote more efficient resource acquisition, leading to higher productivity and higher carbon stocks. The basal area and tree height are both good indicators of biomass and productivity, and hence, carbon storage, at the stand level [18,39,45]. Both of these stand structural features had a positive influence on biomass and total carbon storage in homegardens, according to bivariate relationships and subsequently verified by structural equation models, with the basal area having the strongest effect (Figures 3 and 4a,b,e,g). As a result, in order to maximize and sustain both above- and belowground carbon storage (including soil organic) in homegardens, we must maintain higher diversity and stand tree height, where tree diversity increases stand basal area and improves total carbon storage in homegardens. Recent studies in tropical forests have shown that that both selection and complementarity directed the biomass carbon stock together [13]. Potential work should consider different components of diversity, e.g., species richness, functional trait diversity and functional trait composition, to investigate the mechanism between diversity and carbon relationship in homegardens.

#### *4.5. Implication for Nationally Determined Contributions*

Agricultural activities are among the main sources of carbon emissions in developing countries. Integrated agroforestry systems, however, are considered to be sustainable landuse systems, such as traditional homegardens, which can be a crucial tool for achieving carbon mitigation goals because it has high potential for sequestering atmospheric carbon. Nevertheless, some countries have not included homegardens as possible NDCs, such as Bangladesh, whose homegardens coverage is so rich. Here, we showed that in the Digholia sub-district, each village had high tree coverage (Figure 2) with high diversity (Table 3) and carbon storage (Table 4). In its initial forest reference level, homegardens have not been

listed as forests by the FAO (Food and Agricultural Organization) definition in Bangladesh. Nevertheless, our study shows that both tree coverage (24%) and mean tree canopy height (7.87 m) in homegardens within homesteads satisfied the forest definition by the UNFCCC (United Nations Framework Convention on Climate Change; tree cover: 10–30%; Canopy height: 2–5 m) [46]. Through considering many homegardens as a unit, we can view it as continuous forests (Figure 2). Therefore, homegardens also follow the concept of minimum area of land (0.05–1.0 ha) [42]. Thus, districts with homegardens in Bangladesh satisfy the requirements of REDD+ or the IPCC (Intergovernmental Panel on Climate Change) definition of forests, which has been used in other forest-based mitigation alternatives such as the clean development mechanism (CDM). Nonetheless, since homegardens are not currently included in the REDD+ scheme, we should explicitly enlist homegardens as a choice for prospective NDCs (although a large-scale study is needed in other districts of Bangladesh). Therefore, by setting up a Village Forestry Committee or a group (composed of several homegardens owner), we can manage it like community-based forestry program [18], e.g., social forestry (managed by the Social Forestry Committee with Gender Equity under direct oversight of the Forestry Department of Bangladesh). In this way, it can retain some of the contributions to the Government's climate change adaptation and mitigation measures promised by the Paris Agreement (NDCs; 5% reduction in greenhouse gas emissions by 2030).

#### **5. Conclusions**

In this study, we assessed the tree coverage and direct and indirect effects of tree diversity on carbon storage through stand structure in homestead agroforestry systems (homegardens) in Dighalia, using both remote-sensing and field data. Remote sensingbased tree coverage and consistency of extent of homegardens and field-based canopy height met the forest definition as per IPCC or UNFCCC protocol. High tree diversity with a higher stand structure and carbon storage both above- and belowground indicates that homegardens in Bangladesh play a key role in the conservation of biodiversity and climate change mitigation. We also found that tree diversity promotes carbon storage directly and indirectly via basal area. Therefore, we suggest that for maximizing and sustaining carbon storage, we need to maintain diverse tree species which would lead to a higher stand basal area and, hence, high carbon storage in homegardens. In this way, districts with high homegardens in Bangladesh could contribute to mitigating carbon emissions, and thereby, be a prospective choice (as an agroforestry system) for the NDC.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/f12121661/s1, Table S1: Result of multicollinearity test of different covariates using variance inflation factors analysis. Table S2: Fit statistics of the rejected and accepted structural equation models for predicting different carbon pools in combination of plant diversity and stand structure.

**Author Contributions:** Conceptualization, M.M.R. and G.K.K.; methodology, M.M.R. and G.K.K.; software, M.M.R.; formal analysis, M.M.R. and G.K.K.; investigation, G.K.K. and H.A.; writing original draft preparation, M.M.R. and G.K.K.; writing—review and editing, M.M.R. and M.E.K.; visualization, M.M.R.; supervision, M.E.K. and M.X.; funding acquisition, G.K.K. and M.M.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** Gauranga Kundu provided funding for field data collection and soil organic carbon assessment. This research was also supported by the National Key Research and Development Program of China (2017YFA0604300, 2018YFA0606500).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data analyzed in this manuscript will be provided by the corresponding author upon reasonable request.

**Acknowledgments:** We would like to thank Mahmood Hossain, of Khulna University, Khulna, Bangladesh, for providing laboratory support in Nutrient Dynamic Lab. We want to acknowledge the European Space Agency for providing the image data. We also want to thank Steven Xu from the Johns Hopkins University for improving the English writing of the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

