Landscape-Scale Biodiversity Impacts Analysis of Côte d’Ivoire’s Cocoa Cultivation along Export Supply Chains

Agricultural land use for export commodities leads to significant biodiversity impacts. A spatially detailed assessment of these impacts is crucial for implementing effective mitigation policies. Using cocoa cultivation and exports in Côte d’Ivoire as an example, we present a novel framework that combines earth observations, enhanced landscape-scale biodiversity models, and subnational export supply chain data sets to track the tele-connected potential biodiversity impacts of export groups and importing countries. We found that cocoa cultivation accounts for ∼44% of the biodiversity impacts in Côte d’Ivoire’s cocoa cultivation areas, with >90% attributable to cocoa exports. The top 10 importing countries account for ∼84% of these impacts. Our method offers improved spatial detail compared to the existing approaches, facilitating the identification of biodiversity impact hotspots. Additionally, the biodiversity impacts of agroforestry cocoa are not always lower compared to full-sun cocoa, especially when agroforestry systems are established in regions of high biodiversity importance. Our transferable framework provides a comprehensive understanding of biodiversity footprint and promotes informed decision-making for sustainable agricultural production, processing, and trade. Our framework’s application is currently constrained by the scarcity of detailed supply chain data sets; we underscore the urgent need for improved supply chain transparency to fully unlock the framework’s potential.


Methods S1: Mapping a harmonized land use map for Côte d'Ivoire
Pre-processing of input datasets: We constructed a harmonized land use map of Côte d'Ivoire by integrating 7 raster, vector, and numerical datasets based on remote sensing, and national or regional statistics.To enable the integration of these datasets, we preprocessed the data for consistency.First, all spatially explicit data layers were clipped to the spatial extent of Côte d'Ivoire.The Côte d'Ivoire border vector data are available on the open-source site GADM 1 .We then use the nearest neighbor method to resample the raster or vectors to the spatial extension of the BNETD 2016 map with a spatial resolution of ~10m.

Agriculture allocation:
As our point of interest was the land-use-related biodiversity impacts of cocoa cultivation, we give the highest priority to cocoa grid cells and differentiated cocoa cultivation types.Cocoa cultivation types include full-sun cocoa and agroforestry cocoa.Full-sun cocoa, also known as monoculture cocoa, is typically grown in large, open fields where the trees are exposed to full sunlight.Agroforestry cocoa is grown more diversely and sustainably.In agroforestry systems, cocoa trees are intercropped with other plant species, such as shade trees and fruit trees.We differentiate cocoa cultivation types by integrating cocoa map and BNETD 2016 map.Cocoa grid cells identified as 'Agroforestry management' and 'Degraded agroforestry management' in BNETD 2016 map are allocated to agroforestry cocoa, while other cocoa grid cells are allocated to full-sun cocoa.
We allocated the other agricultural grid cells from BNETD 2016 map.The grid cells labeled as 'Coffeecocoa', 'Rubber tree', 'Coconut palm', and 'Cashew tree' were allocated to 'Low-input agriculture'.
The grid cells classified as 'Agroforestry management' and 'Degraded agroforestry management' were allocated to 'Agroforestry'.The agriculture grid cells in Côte d'Ivoire were not classified as 'Highinput agriculture' because we used the GLOBIO4 2 methodology to calculate land use intensity.Based on the results, all agriculture grid cells in the country are considered low input.The results is also consistent with existing research's assumptions 3 .

Pasture allocation:
We took the average value of the FAOSTAT land use statistics for 'Land under temp.meadows and pastures' and 'Land under perm.meadows and pastures from 2015 to 2020 as the total pasture area.
We allocated pasture grid cells according to their suitability values and the total pasture area from the grid cells with the highest suitability until the total area was allocated.Tropical Livestock Units (TLUs) are a way of standardizing and comparing different types of livestock based on their weight and size in tropical regions.We first computed each grid cell's total TLUs based on the density of ruminants (cattle, goat, sheep) in the Gridded Livestock of the World database (GLW v4) 4 datasets.The conversion factors we used are one goat or sheep = 0.1 TLUs and one cattle = 0.6 TLUs 5 .We specified that pasture could only be allocated to grid cells categorized as 'Wooded savannah', 'Shrub formation', and 'Herbaceous formation' from BNETD 2016 map.We utilized Potential Natural Vegetation (PNV) map to differentiate between 'Man-made pastures' and 'Livestock grazing'.Pasture grid cells with herbaceous/shrub vegetation PNV are categorized as 'Livestock grazing'.However, since all allocated pasture grid cells had tree PNV, we classified all pasture grid cells as 'Man-made pastures'.

Forest allocation:
We allocated the remaining forest grid cells in BNETD 2016 map to different types of forests.We used the ~100m spatial resolution Global Forest Management Data (GFMD) 6 map to obtain the spatial distribution of primary and plantation forests.GFMD map divides the other forest grid cells into naturally regenerating forests with signs of human activities, oil palm plantations, and agroforestry.
We did not use these three forest types for allocation because agroforestry and oil palm plantations have been considered in agricultural land use.We used the Forest Landscape Integrity Index (FLII) 7 for further allocation.FLII is a globally consistent, continuous index of forest condition as determined by the degree of anthropogenic modification by integrating data on observed and inferred human pressures and the index of lost connectivity.We allocated the remaining forest grid cells with FLII larger than 0.6 as 'Lightly used natural forest' and the FLII smaller than 0.6 as 'Secondary forest', corresponding to medium integrity and low integrity in the original paper.The remaining unallocated grid cells are allocated to 'primary vegetation', which is mostly present in protected areas.

Methods S2: GLOBIO-InVEST modeling
The GLOBIO model provides a biodiversity intactness index based on mean species abundance (MSA), which is the average response of populations of different species to various stressors such as land use changes, fragmentation, and infrastructure 8 , MSA ranges from 0 to 1, where 1 means the species composition is completely intact, and 0 means that all original species are extirpated (locally extinct).
The GLOBIO-InVEST model extends the GLOBIO3 method to downscale their global approach to the landscape level to identify finer ecological responses that may include nonlinearities 9 .The main differences between GLOBIO-InVEST and GLOBIO3 are that GLOBIO-InVEST uses a more sophisticated approach to quantifying fragmentation than applying overall averages of patch size in different habitats.GLOBIO-InVEST replaces the standard GLOBIO3 model with fragmentation analysis using a fragmented forest quality index (FFQI).The FFQI is calculated by considering how many of a forest's neighboring cells are also forested.GLOBIO-InVEST MSA simulation is based on the same cause-effect philosophy as GLOBIO3.The details of the cause-effect parameters we used in this study are in Table S3.The parameters are based on a broad literature review and suggested methodologies for inputting and processing the spatial data required.During the simulation, suppose the land use type of grid cell i is lu, the proximity of infrastructure is I, and FFQI is F. The MSA for grid cell  is calculated by multiplicative as follows： We assumed that the cause-effect parameters satisfy a normal distribution to obtain the upper and lower bounds for the MSA simulations.We used the mean parameters plus or minus two times its standard deviation to determine the MSA distribution with 95% confidence intervals.The range of the MSA cause-effect parameters are limited from 0 to 1.

Methods S3: Comparative Analysis
We compared the results of the BIM method with the PDF method to illustrate the perspective brought by our landscape-scale assessment approach.PDF accounts for a fraction of species richness potentially lost due to an environmental mechanism.We employed the UNEP life cycle initiative recommended CFs for PDF calculation.The CFs were generated using the countryside species−area relationship (c-SAR) and vulnerability scores to estimate PDF per unit area of land occupation/transformation in 804 terrestrial ecoregions across five taxa and six land use types 10 .We only considered the land occupation of cocoa cultivation, and the CFs we used for the different ecoregions are shown in Table S4, and the spatial distribution of ecoregions in Côte d'Ivoire is shown in Figure S7.We first rasterized the ecoregion vector files to the spatial extension of the satellite remotely sensed cocoa map.We then allocated CFs based on the ecoregion corresponding to each cocoa grid cell.Finally, we multiplied the CFs of each cocoa grid cell with the grid cell area to obtain the global PDF of each cocoa grid cell.We calculated the PDF per ton of cocoa produced and linked PDF results to the export supply chains in the same way as in sections 2.2.3 and 2.3.to man-made pastures, a figure that maintains consistency with the FAOSTAT data since we use FAOSTAT data as the total area reference for pasture.
There is a clear discrepancy between the GLOBIO4 and HILDA+ land use maps in southern Côte d'Ivoire.GLOBIO4 typically classifies most of the grid cells in this region as agricultural land, while HILDA+ classifies them as forest.Therefore, instead of grouping tree grid cells into a single category, our map effectively distinguishes perennial tree crops and agroforestry from tree cover by integrating detailed local land cover maps and allocation rules.Comparatively, the total agricultural and forest land allocation of our map lies between the results of GLOBIO4 and HILDA+ and is closer to the FAOSTAT reporting data (Figure S9).
Through visual inspection using Sentinel-2 cloudless 2020 and the high-resolution land cover map (WorldCover 2020), we further evaluated these land use maps in typical landscapes (Figure S8).The visual inspection shows that our harmonized land use map provides more detailed spatial granularity and landscape information than both GLOBIO4 and HILDA+.
The discrepancies between existing land use maps and the precision offered by our approach underscore the complexity and dynamic nature of land use classification.Moreover, our approach demonstrates the importance of integrating multiple data sources and validation techniques to produce more reliable land use maps.It underscores the utility of detailed support layers and locally sourced land cover data in refining land use datasets.

Results S2: Comparison of cocoa cultivation maps
The total remotely sensed cocoa cultivation area was 4,450 kha, ~0.5% less than the FAOSTAT figure.
Agroforestry cocoa accounted for ~30% of the total cocoa cultivation area.The cocoa cultivation area aligns well with climatically suitable growing regions 13 and dominates land use in these regions.In the cocoa climatically suitable growing regions, over 26% of the land is used for cocoa cultivation (Figure 2b and S10).In Bas-Sassandra, cocoa cultivation even accounts for ~42% of its total land area.
We compared the remotely sensed cocoa cultivation map used in this study with the Spatial Production Allocation Model (SPAM) 14 and another remotely sensed cocoa map (referred as UoW) 15 .To facilitate comparison, we resampled both remotely sensed cocoa distribution maps to the spatial extent of SPAM.
First, there is a significant difference in the cocoa cultivation maps.The difference of the cocoa cultivation area can vary by four times in the same region (e.g., Bas-Sassandra).Second, the cocoa distribution in both remotely sensed maps was more concentrated in the southwest and southeast regions of Côte d'Ivoire, which corresponds well with the climatically suitable growing regions 13 , while the SPAM cocoa map showed a more dispersed distribution (Figure S11).This dispersed distribution caused its area in cocoa climatically suitable growing regions being much lower than the remotely sensed maps.This discrepancy will prevent us from correctly identifying the biodiversity impact hotspots of cocoa cultivation.
We further verified the cocoa cultivation map by assessing the consistency between cocoa production estimated based on cocoa map and cocoa production data from Trase supply chain.For the convenience of comparison, we first aggregated the cocoa production into each department.The production estimates from the remotely sensed cocoa map showed lower values than Trase records export volume in only 5 departments.The total negative deviation is approximately 12,303 tons.On the other hand, SPAM produced similar results in 14 departments, showing a total negative deviation about 18 times higher than that based on remotely sensed cocoa maps.We also tested the correlation of cocoa production of the spatially explicit part of the Trase dataset with production estimates from the remotely sensed cocoa map (Figure S12a) and with SPAM (Figure S12b).We found that the relation with remotely sensed cocoa production estimates (R = 0.81) is higher than with SPAM (R = 0.48).The above comparisons suggest that the remotely sensed cocoa map provides a more accurate representation of cocoa cultivation in Côte d'Ivoire and is more consistent with Trase supply chain data.

Results S3: Uncertainty analysis
In Côte d'Ivoire, we found that the average relative errors of the biodiversity impacts of cocoa cultivation from using different land use maps (specifically HILDA+ and GLOBIO4) and different cause-effect parameters varied significantly, ranging from -62% to 101%.This error is much smaller than the PDF method (the average relative error of the PDF method is more than 400%).We performed further checks, and the relative error does not have a clear impact on the identification of biodiversity impact hotspots and our other results.We have detailed the upper and lower bounds/ relative errors of the biodiversity impacts of cocoa cultivation in each cocoa production department in Table S8.The upper and lower bounds of the biodiversity impacts per ton of cocoa produced in each cocoa production department are detailed in Table S9.In addition, the upper and lower bounds of the tele-connected biodiversity impacts of cocoa imports in different importing countries are listed in Table S13.

Figure S2 .
Figure S2.(a) The spatial distribution of the aggregated rarity-weighted richness.(b) The spatial distribution of the rarity index (RI).(c) Histogram of the frequency distribution of the aggregated rarity-weighted richness.(d) Histogram of the frequency distribution of RI.

Figure S3 .
Figure S3.Diagrammatic representation of the biodiversity indicators utilized in this research.

Figure S4 .
Figure S4.The BNETD 2016 land cover map used to create the harmonized land use map in this study.

Figure S5 .
Figure S5.The GLOBIO 4 land use map (2015, 300m).Forest or shrub grid cells that are not allocated to human activity land use types are integrated as 'Unallocated vegetations'.

Figure S9 .
Figure S9.Area Statistics of land use in Côte d'Ivoire, with the harmonized land use map for this study represented by stacked bars, and HILDA+ and GLOBIO4 by a single bar, the data reported by FAOSTAT are represented by black horizontal lines and numbers.

Figure S10 .
Figure S10.Area statistics of different cocoa maps across climatically suitable growing regions, with black horizontal lines indicating the percentage of the region's total land area under cocoa cultivation.

Figure S11 .
Figure S11.Comparisons of the cocoa cultivation maps.(a) remotely sensed cocoa map 16 used in this study, (b) remotely sensed UoW cocoa map 15 and (c) SPAM cocoa harvested area map 14 .The cocoa maps were all resampled to the SPAM grid cells for easier comparison.d-f shows the relative differences between the different cocoa maps.

Figure S12 .
Figure S12.Correlations between the cocoa production estimates with Trase export volume at each department.(a) Correlation graph between cocoa production estimates based on remotely sensed (RS) cocoa map and Trase export volumes; (b) Correlation graph between SPAM cocoa production estimates and Trase export volumes.

Figure S13 .
Figure S13.(a) the spatial distribution of the PDF of cocoa cultivation.(b) the PDF per ton of cocoa produced in each department.

Figure S14 .
Figure S14.Comparative analysis of biodiversity impacts and biodiversity impacts per ton of cocoa produced in the top 10 cocoa production departments.Biodiversity impacts, or biodiversity impacts per ton of cocoa produced, is presented as the relative percentage of the total biodiversity impacts due to cocoa cultivation in Côte d'Ivoire for easier comparison.The bar chart is used to represent the relative biodiversity impacts of each department.The circle represents the relative BIM per ton of cocoa produced by each department, and the triangle represents the relative PDF per ton of cocoa produced by each department.

Figure S15 .
Figure S15.Panels a-c display box-and-whisker plots, illustrating the biodiversity indicators for each department with production volumes surpassing 1000 tons.The biodiversity indicators are as follows: (a) production weighted average MSA loss; (d) production-weighted average RI; (e) BFcocoa.The box-and-whisker plots' terminal ends denote 95% confidence intervals, while the filled circles symbolize the median of the biodiversity indicator.Solid black triangles denote the mean biodiversity indicator.The dashed circles represent the indicators of each department.

Table S1 .
Definition of land use types in this study.Secondary forest A forest that has re-grown after a timber harvest or clearing for agriculture, until a long enough period has passed so that the effects of the disturbance are no longer evident.
6Man-made pastures The pastures are created by removing natural vegetation and then sowing grasses or other forage crops that are suitable for animal grazing.7AgroforestryAgroforestry is a land use type that integrates the cultivation of trees or shrubs with crops or livestock on the same piece of land.The 'agroforestry' in this study was limited to land use for cocoa cultivation under the canopy shade above 15m.8Low-inputagricultureLow-inputagriculture is a land use type of farming that relies on minimal external inputs, such as fertilizers, pesticides, and energy, and instead emphasizes the use of natural resources and ecological processes to maintain soil fertility and manage pests and diseases.10Built-upareasThepresence of buildings (roofed structures).Includes paved surfaces (roads, parking lots), commercial and industrial sites (ports, landfills, quarries, runways), and urban green spaces (parks, gardens).* In the biodiversity intactness modeling, 'Primary forest' and 'Primary vegetation' were combined into one land use type.

Table S2 .
The input datasets used to construct the land use map of Côte d'Ivoire.

Table S3 .
The cause-effect parameters used in modeling MSA.

Table S4 .
The characterization factors (CFs) used for potential disappeared fraction (PDF) calculation.

Table S5 .
Land use allocation results and comparison with other land use datasets.

Table S6 .
Departmental statistics of cocoa cultivation area in Côte d'Ivoire.
* Only includes departments with an area greater than 1Kha.

Table S7 .
Departmental statistics of cocoa production in Côte d'Ivoire.
* Only departments with cocoa production exceeding 1000 tons are included.

Table S8 .
Departmental statistics of biodiversity impacts of cocoa cultivation in Côte d'Ivoire.

Table S9 .
Departmental statistics of land-use-related biodiversity impacts per ton of cocoa produced in Côte d'Ivoire.

Table S10 .
Statistics of land-use-related biodiversity impacts and biodiversity impacts per ton of cocoa imported for importing countries.

Table S11 .
Departmental statistics of biodiversity impact indicators in Côte d'Ivoire.The biodiversity impact indicators included as follows: production weighted average MSA loss, production-weighted average RI, and BFcocoa.
* Only departments with cocoa production exceeding 1000 tons are included.

Table S12 .
Statistics on the biodiversity impacts of Dutch cocoa imports from Côte d'Ivoire for each cocoa production departments in absolute and relative perspectives (only includes the share that can be traced back to specific production departments).