2.1. Study area
We carried out this study in farms and municipalities that produced C. arabica, distributed in the Northwest, North-Central, and “Pioneiro” North regions, located on the northern edge of the state of Paraná, Southern Brazil (Figure 1). The native vegetation cover represents only 11% of land use (MapBiomas, 2021). The region currently corresponds to a matrix of intense agricultural activity, with emphasis on the production of corn (Zea mays) and soybean (Glycine max) and represents the southern limit of coffee plantations in the country. The sampling units are located within the Atlantic Forest biome in a region characterized by the humid subtropical (Cfa) and oceanic temperate (Cfb) climate domain, with a phytogeographic domain of mixed ombrophilous and semideciduous seasonal forests, and an average elevation of 653 m (±122). In this region, C. arabica is the main species cultivated for coffee beans production.
2. 2. Data survey
We obtained data on coffee production per farm from interviews with 30 coffee producers in the northern region of the state of Paraná, southern Brazil, who had their contacts provided by local cooperatives. We selected data from farms with more than 1 ha Arabica coffee-growing areas in 2018, totaling 25 sampling units. We then estimated the production data per area, using the ratio of kilograms per hectare, as follows. In the interviews, the mass production of farms was obtained in bags of green beans through production records, and later converted into kilograms. We estimated the production area considering only the productive stands, according to the owners, and using the QGIS v3.16 Hannover program (QGIS Development Team, 2020). We calculated the area by crossing remote sensing images taken with the Landsat8 satellite, from July to December 2017 (INPE, 2018), with the maps of coffee-growing areas for the state of Paraná in the year 2017 (CONAB, 2018b). We obtained production data by municipality from the most recent agricultural census (IBGE, 2017). We used information from the 11 municipalities in which the 25 farms were distributed. To increment our data, we added another 19 municipalities drawn, totaling a sample of 30 municipalities. We included only municipalities with more than 50 coffee-growing farms. For the sampled municipalities, we extracted the variables of production (tons) and area of coffee harvested (hectares), and thus estimated the production per area (kg ha^-1).
2.3. Land use classification
We used the land use classification from the 6.0 collection of the MapBiomas Brazil Project (2021) for the state of Paraná in the year 2017, the year that corresponds to the flowering of the 2017/2018 harvest and the last agricultural census carried out in Brazil. The land use classification of this collection was generated through the Random Forest algorithm from remote sensing images produced by the Landsat8 satellite, which has a pixel definition of 30 x 30 m (MapBiomas 2021). We considered all the 15 classes described in MapBiomas occurring at the samples to calculate landscape metrics. We performed a manual correction of landscape classes to eliminate possible discrepancies with more detailed mappings and satellite image observations. For this end, we used the mapping of forest remnants from SOS Mata Atlântica (2019) and remote sensing images from the Landsat8 satellite from July to December 2017 (INPE, 2018). We processed the images and classified the landscape using the “Semi-Automatic Classification Plugin” (Congedo, 2016) in the QGIS v3.16 Hannover program (QGIS Development Team, 2020).
2.4. Landscape metrics
To assess the effect of landscape composition and configuration on coffee production, we calculated six landscape metrics for both landscape and municipal levels. We calculated only the native forest class for class-level metrics, due to the importance of this class as an ecosystem services supply area for coffee crops (Ricketts et al., 2004; Hohlenwerger et al., 2022). To measure the density of forest fragments, we used the density of native forest fragments, which corresponds to the number of patches of native forest in an area of 100 ha. As a measure of forest isolation, we used the average Euclidean distance from the nearest neighbor of the class. To measure the proximity of forest patches to farms, we used the average proximity index of forest patches to the central point of the buffer (centroid of the farm area). We measured the connectance index of native forest patches within the specified threshold distance of 2,000 m. We used this distance, considering the average flight radius of bees, the main pollinators of coffee (Saturni et al., 2016). In this way, pairs of forest fragments with a linear distance < 2000 m from the centroid of the patches were considered connected. To measure the diversity of the landscape we used the Shannon equitability index, which weights the occupation of a class by the total of classes, thus measuring the relative distribution of coverage of each land-use class. Finally, to assess the level of anthropogenic intensity in the landscape, we used the landscape intensity index, which results from the proportion between the anthropogenic landscape (classes of agriculture, livestock, and urban infrastructure) and the natural landscape (classes of natural areas) (Hipólito et al. al., 2018). To calculate this metric, we did not include the classes corresponding to water bodies.
We calculated the landscape metrics for radii from 1 to 5 km from the centroid of the farms at 1 km intervals. We defined the radii based on the displacement and foraging radius of species that provide essential ecosystem services for coffee growing, such as insectivorous birds and bees (Buler et al., 2007; Saturni et al., 2016). We calculated all the metrics in the FRAGSTATS v4.2.1 program (McGarigal et al., 2012) except for the landscape intensity index calculated in the R Studio environment (R Core Team, 2021). We then tested the autocorrelation among the predictor variables (landscape metrics) calculated in the 5 km radius and we excluded them when the correlation index was greater than 0.70. Forest isolation and forest proximity index (Appendix A Figure S2) were excluded from the analysis since their correlation with the response variable (R²) was weaker than their competitors (Table 1). Finally, we tested the scales of effect of the variables from linear models (LM) performed with the “lm” function (package “stats”), where the response variable was related to each predictor variable for each radius of influence. The variables were modeled only on the radius with the greatest effect (Miguet et al., 2016, Table 1).
Table 1 Correlation coefficient (R²) between landscape metrics and coffee production in farms in Southern Brazil, Paraná State, Brazil, for the different radii of analysis. The highest R² value obtained for each variable is highlighted.
|
Effect radius
|
Variable
|
1Km
|
2Km
|
3Km
|
4Km
|
5Km
|
Forest cover percentage
|
0.132
|
0.139
|
0.144
|
0.119
|
0.104
|
Forest patches density
|
0.122
|
0.273
|
0.235
|
0.170
|
0.188
|
Forest proximity index
|
0.049
|
0.019
|
0.006
|
0.012
|
0.017
|
Forest isolation
|
0.226
|
0.357
|
0.323
|
0.396
|
0.392
|
Forest patches connectance
|
0.173
|
0.009
|
0.108
|
0.143
|
0.046
|
Landscape diversity
|
0.203
|
0.231
|
0.213
|
0.213
|
0.198
|
Landscape intensity index
|
0.388
|
0.400
|
0.392
|
0.357
|
0.354
|
2.5. Statistical analysis
To ensure our data on coffee production was robust face to seasonal and regional variations and represented the Brazilian production, we compared the average production per area of farms and municipalities to the average Brazilian production (CONAB 2018a) based on a t-test for one sample with the “t.test” function (package “stats”). We established a standardized value of coffee production for each farm from the difference between the average Brazilian production and the production of either farms or municipalities. We used the standardized value as a reference to verify the relationship between production variation and the surrounding landscape. To show the variation of farm/municipal production according to relevants landscape metric, we plotted the standardized values on a graph intercepted by a horizontal line corresponding to the upper quartile of the variation in coffee production (from 417,680 to 872,721 kg ha^-1) to indicate the most productive farms.
We verified the spatial autocorrelation of the response variable through a Mantel test with the “mantel” function (vegan package), considering farm coordinates and the dependent variable (Appendix A Figure S3). To verify how coffee production (kg ha^-1) varies as a function of the density, isolation, connectivity and proximity of native forest patches, the diversity and level of landscape intensity index, we carried out a selection of univariate linear models generated from the “glm” function (package “stats”) (Johnson et al. 2014). The models were structured with a Gaussian distribution and identity link function and validated through the criteria of normality, homogeneity and outliers of residuals through the “TestResiduals” function (package “DHARMa”). Once the best univariate model was detected, we performed a new selection of models, this time from multivariate models. The models were structured by pairwise combining the predictor of the best univariate model in interaction with the other predictors (Johnson et al. 2014). Model selection was performed using the second-order Akaike information criterion (AICc) through the “model.sel” function (package “MuMIm”), where the best model was determined by the smallest AICc value. We considered plausible models with ΔAICc ≤ 2, AICc higher than the model with no effect and weight ≥ 2.0 (Cavanaugh and Neath, 2018). We performed statistical analyzes using the RStudio v.1.2.1 statistical programming environment (R Core Team, 2021). The analysis scripts and data used in this study are available in the GitHub repository (directory “DataCoffePR”).
2.6. Ethical issues
In order to carry out the interviews with the coffee farmers, we submitted this research project d to the Humans Research Ethics Committee (CEP/SD) of the Federal University of Paraná and approved under CAAE: 13366919.2.0000.0102.