Higher bee abundance, but not pest abundance, in landscapes with more agriculture on a late-flowering legume crop in tropical smallholder farms

Study area and field selection

We conducted our study from May to August 2019 in Mzimba district, Northern Malawi. We selected ten already existing pigeon pea fields. We were granted verbal permission for conducting the research on each of the farmer’s private fields. Their names are: Isobel Lubanda, Adams Tembo, Mercillina Tembo, Ireen Mhoni, Simon Chitaya, Jacob Mvula,Jane Salanda, Lyna Njunga, Goodson Moyo, Moles Thupa. The farmers are not represented by a company or a farming cooperation, but were in contact with the authors through the SFHC (Soils, Food and Healthy Communities) organization. We have no form of written permission for the conduction of the research. In all the ten fields, the pigeon pea crop had been planted at the onset of rains in December 2018, and were initially intercropped with groundnut (Fabaceae: Arachis hypogaea L.). By the time we began data collection, all the groundnut had already been harvested from all the fields. All pigeon pea fields selected were planted with a local medium-maturing variety. The peak of bloom of this pigeon pea variety is in May in our system. The duration of the flowering period can depend on the climatic conditions, but in our region, the bloom lasted about 4 weeks.

Malawi is located in the seasonal tropics and experiences a marked peak in rainfall from December until the end of February. In the months when we performed our experiment, there was no rainfall, as is typical during this time of year (Mungai et al., 2016). The pigeon pea in our study region is a rain-fed crop and is not irrigated or watered in any way, especially as pigeon pea is considered drought-resistant (Odeny, 2007). All field activities, including land preparation and weeding, were managed traditionally by hand. Farmers did not apply any chemicals such as herbicides and pesticides on their fields.

Fields ranged from 166 m ² to 577 m² in size, with mean field size being 332 m². This is representative of field sizes of such a crop in the study region, where the average smallholder total farm size ranges from only 0.5 to 1.4 hectares (FAO, 2018). Field size did not correlate significantly with the proportion of semi-natural area (F_2,6 = 2.08, R² = 0.21, p = 0.683) nor with the proportion of agricultural area (F_2,6 = 2.08, R² = 0.21, p = 0.088) in the 1 km radius surrounding our fields. Field margins were vegetated with non-flowering weeds, grass or shrubland. As it was the dry season during data collection, there were no flower margins on the fields. The surrounding agricultural fields where mostly empty, as the main staples in Malawi, such as maize, was already harvested by this time in the season. Surrounding semi-natural habitat was mainly composed of shrubland and forest. Generally, these are not actively managed but may to some extent be exposed to exploitation by people due to economic activities such as collection of firewood and grazing of livestock.

We aimed to choose sites which were at a distance of at least 2 km from each other. However, one site was found to have too large an overlap with two others within a 1 km radius, with the center of this field being 883 and 885 m away from the center of the nearest and second-nearest site, respectively. Therefore, this site was subsequently dropped from any landscape analyses (Fig. 1). The remaining fields were located within two non-correlating gradients of semi-natural habitat (ranging from 2% to 32%), and agricultural area (ranging from 25% to 75%) within a 1 km radius surrounding the fields (F_1,7 = 0.56, R² = 0.07, p = 0.480). The 1 km radius was chosen, because we wanted our sites to be independent from each other and prevent spatial autocorrelation. Additionally, since bees are central place foragers, and their foraging ranges are limited, most bees would be sensitive to landscape differences within the 1 km radius (Steffan-Dewenter et al., 2002; Zurbuchen et al., 2010). Other habitats in our study area included some built-up areas (such as buildings and roads) and bare rock (mostly on hilltops). Although honeybees are native to the area (Requier et al., 2019), we found no honeybee hives in any of the fields across our study area. Moreover, none of the farmers we worked with kept honeybees on any of their fields. In our study area, honeybees are rarely actively placed in fields by farmers, but rather encouraged to nest nearby by placing traditional beehives near fields where they may be passively colonized by a honeybee colony (Appendix 2) (Requier et al., 2019). To our knowledge, there were no such hives placed near any of our study fields.

Landscape analysis

For the land use and land cover classification, we acquired three cloud-free Sentinel-2 satellite images from 2019 from the Copernicus Open Access Hub (https://scihub.copernicus.eu/dhus/#/home). One image was taken on November 8th 2019, and two are from November 15th 2019, which is in the late dry season in Malawi. Though this is not the same time as our field study period, the images still show the general land use and land cover in our study area and we assume this has not changed significantly between May and November.

The methodology for land cover analysis involved images pre-processing, supervised maximum likelihood classification, post-classification, and accuracy assessment (Appendix 3). First, we pre-processed the Sentinel-2 images, which consisted of the atmospheric correction, image resampling, layer stacking, seamless mosaic, and image subsetting. Then, we conducted Maximum Likelihood classification using training samples generated from fieldwork and Google Earth. The classification includes six classes, which are: (1) bare land/road; (2) shrubland; (3) agricultural land; (4) water/riverbed; (5) settlement; (6) trees/forest. However, in this study, we only used classification (2) shrubland and (6) trees/forest together which we defined as semi-natural habitat (SNH) and (3) agricultural area. Finally, we performed post-classification and accuracy assessment. The overall accuracy of the classification is 85.1%, with a Kappa Coefficient of 0.817. We used tabulate area tools in ArcGIS to get the area and proportion of the types of land use and land cover for all buffer zones of each site (ArcGIS Pro 2.6; Esri, Redlands, CA, USA).

Bee and beetle collection

To assess bee and blister beetle densities, we walked straight 15 m transects for 10 min across the most flower rich area of the field on three separate dates, resulting in 30-minute total sampling time per field. One round of transects at each field was finished before a new round was started, resulting in a pause of about a week between individual sampling dates at each field. All flower-visiting bees one-meter left and right of the investigator were captured with an insect net for subsequent identification. All blister beetles observed one meter left and right of the investigator were counted, taking care not to count the same individual twice. To obtain a proxy for planting density, we counted the number of pigeon pea shrubs across each transect. Planting density across transects did not correlate with the proportion of agricultural area (p = 0.786) or semi-natural habitat (p = 0.338) (F_2,7 = 0.53, R² = 0.13, p = 0.610). Transects were walked between the 8th and 29th of May 2019, between 8:00 and 16:00, when weather conditions were not windy (Beaufort scale <  3) or too cloudy (<80% cloud cover), and thus deemed favorable for bee foraging. It did not rain, and temperatures ranged from 18 °C to 25 °C throughout our sampling period. Visitation order on a specific day was randomized so that each field was visited during a different time of the day during consecutive transects.

Bee identification

Captured bees were identified to genus or subgenus level with the guides from Michener (2007) and Connal, Kuhlmann & Pauly (2010) and grouped by (morpho)species. Captured bees are stored at the Biocentre, University of Würzburg, Germany.

Flower exclusion and fruit set data

At each site, we marked 15 pigeon pea plants in a continuous line from the edge of the field inwards. On each plant, we tagged one cluster of flowers as the open control. This cluster was accessible to all visitors, both pollinators and herbivores. On the same plant, we then paired this tagged cluster with another cluster of flowers to which all visitors (pollinators as well as herbivores) were excluded using a 9 × 12 cm organza bag. The number of flowers in the tagged and bagged clusters where counted. The bags had a mesh size of 0.6 mm—small enough to exclude any insect. Although exclusion of pollinators and pests in different treatments would have been ideal, this could not be done, since the beetles feed on the flowers, during the same time that pollinators are visiting them. Our hypothesis was that if the fields experience high pest pressure, bagged clusters will perform better, as they are protected from herbivory. On the other hand, we assumed that in fields where there is a large amount of ambient pollination and low flower damage by herbivores, the open clusters would outperform the bagged ones. In fields with little ambient pollination, or where the benefits of pollination are cancelled out by pest damage, open and bagged flowers would perform similarly. Plants were tagged and bagged upon the first visit to the field before the flowers had opened and we removed the bags when all the flowers in the cluster had finished blooming, which took approximately two weeks. After removing the bags, the pigeon pea pods were left to mature in the field.

Fruit set data was collected from the 3rd of July to the 2nd of August 2019. To assess fruit set (the proportion of flowers turning into pods) as a measure of pollinator effectiveness, we counted the number of flowers that were originally present on the tagged clusters, and then counted the number of pods formed in the same clusters. The number of pods formed divided by the number of flowers was taken as a measure of fruit set per cluster. In one field, damage by cattle grazing destroyed the tagged plants and we were unable to collect data on fruit set there.

Blister beetle damage assessment

To get a measure of the proportion of flowers damaged by blister beetle herbivory/florivory, we assessed flower damage three times at eight of the sites and twice at two of the sites. We used the open cluster of the 15 pigeon pea shrubs we tagged in each field for this. We counted the number of flowers per cluster and the number of these flowers that showed signs of chewing herbivory typical of blister beetles. With this data we calculated the proportion of flowers damaged by blister beetles.

Data analysis

To test whether landscape composition affected bee and blister beetle abundance, we summed the number of individuals across all three transects. For bee richness, we used the cumulative bee richness across dates per field. We first tested if bee abundance, bee richness and blister beetle abundance were independent of planting density across transects or field size (Appendix 4). We then tested how the proportion of semi-natural habitat and agricultural area within the 1 km radius affected bee and blister beetle abundance using a linear regression. To test whether landscape composition affected bee richness, we used the bee richness at each site and again tested this against the landscape variables using a linear regression. To test to what extent our patterns where driven by the presence of honeybees (Apis mellifera L.), the most abundant pollinator in our system, we tested bee abundance against the two landscape variables including and excluding honeybees from the analysis.

To test if landscape variables affected blister beetle damage in our fields, we calculated for each of the 15 plants the mean proportion of flowers damaged by herbivory across the flowering season from the three dates. Since our data were zero-inflated (no flowers damaged), we used a negative-binomial mixed model using the ‘glmer.nb’ call from the package ‘lme4’ (Bates et al., 2019). We tested the mean proportion of flowers damaged against the proportion of semi-natural habitat and agricultural area. Since we had repeated measures within fields, we used field as a random factor in this model.

To test whether landscapes affected the differences in fruit set between bagged and open clusters, we calculated the proportion of flowers that set fruit for each cluster. Then, we subtracted the proportion of fruit set of the bagged cluster from that of the open cluster. Again, using the package ‘lme4’ (Bates et al., 2019), we calculated mixed effects models testing the difference in fruit set against the two landscape variables, using field as a random factor to account for repeated measures per field. In this analysis, we had to exclude 31 out of 120 plants due tampering or missing tags.

Finally, to calculate the effect of bee visitation and beetle damage on fruit set difference between the bagged and open clusters, we used mixed models. To do this, we took the total number of bees recorded at each site, and divided this by the total number of pigeon pea shrubs across our transects. This gives us bee density per crop plant, which we used as a proxy for bee visitor density per pigeon pea shrub. We then calculated the effect of bee density, bee richness and the proportion of flowers damaged per plant on the difference in fruit set between the paired clusters, using field as a random factor to account for nestedness. We chose to use flower damage, rather than beetle abundance or density in this model because we deemed it a more concrete representation of the pest pressure the plants experienced, though blister beetle abundance and the proportion of damage was correlated (F_1,418 = 4.88, R² = 0.01, p = 0.028). In this analysis, we had to exclude 34 out of 135 pigeon pea shrubs due to tampering or missing tags.

All models were tested for and met the assumptions of distributions, normality (of residuals) and heteroscedasticity. All statistical analyses were performed in R version 4.0.1 (R Core Team, 2020).

Landscape effects on bee abundance, bee richness, and blister beetle abundance

In total, we observed 84 bees of 13 species (Appendix 5) and 127 blister beetles across the five hours of transects during our study period. The proportion of semi-natural habitat within a 1 km radius of the fields did not affect bee abundance (F_2,6 = 5.53, R² = 0.65, p = 0.775) (Fig. 2A) or richness (F_2,6 = 1.38, R² = 0.32, p = 0.203) (Fig. 2C), and neither did it affect blister beetle abundance (F_2,6 = 1.58, R² = 0.35, p = 0.538) (Fig. 2E). The proportion of agricultural area positively affected bee abundance (F_2,6 = 5.53, R² = 0.65, p = 0.0209) (Fig. 2B), though this pattern was primarily driven by higher honeybee densities at high-agricultural area sites, as solitary bees alone did not respond significantly to landscape factors. The pattern was additive, as honeybee densities alone did also not show significant patterns, and it was just the analysis with honeybees and solitary bees together that showed a result (Appendix 6). However, agricultural area did not affect bee richness (F_2,6 = 1.38, R² = 0.32, p = 0.683) (Fig. 2D) nor blister beetle abundance in the fields (F_2,6 = 1.58, R² = 0.35, p = 0.127) (Fig. 2F) (Table 1).

Landscape effects on blister beetle damage and fruit set difference

The proportion of flower damage ranged from zero to 0.36, with a mean of 0.06. There was no effect of the proportion of semi-natural habitat nor agricultural area on the proportion of flowers damaged by blister beetles on the tagged open clusters (Appendix 7). The number of open clusters that a higher proportion of damage than 0.05 varied from 2 to 12 clusters per site, but this did not correlate to either proportion of semi-natural habitat (F_2,6 = 0.01, R² = 0.001, p = 0.923) or agricultural area (F_2,6 = 0.01, R² = 0.001, p = 0.998). The proportion of fruit set on open clusters ranged from 0 (none of the flowers set fruit) to 1 (all flowers set fruit) with a mean fruit set proportion of 0.37. The proportion of fruit set on bagged clusters also ranged from 0 to 1 with a mean fruit set proportion of 0.26. Fruit set difference (open-bagged) ranged from −1 to 1 and had a mean of 0.11. There was no effect of the landscape variables on the fruit set difference between the open and bagged flower clusters (Table 2) (Appendix 7).

Effects of bees and blister beetle damage on fruit set difference

We observed no effects of bee density or bee richness nor of blister beetle damage on the proportional difference in fruit set (Table 3) (Appendix 8).