Study design
This study was performed in the Cantareira-Mantiqueira Corridor (CCM), Atlantic Forest biome located in São Paulo state, Brazil (Fig. 1). The study region is since 2020 a Long-Term Ecological Research (LTER) site, an initiative of the Brazilian federal government that is a network of reference sites for scientific research on the topic of Ecosystem Ecology (CNPq). LTER CCM is a high priority conservation area since it serves as a biodiversity corridor between Cantareira and Mantiqueira hills, two highly forested regions; it also is part of the Atlantic Forest Biosphere Reserve declared by UNESCO in 1993 and, houses the Cantareira system responsible for the water supply of São Paulo city. Remaining forest fragments at LTER CCM are mainly composed of secondary forests that regenerated naturally from coffee and tea plantations abandoned in the late nineteenth century, as well as abandoned pastures in the late twentieth century. The main land uses nowadays are pasture for beef cattle, Eucalyptus plantations, small farms with crop plantations and, urban areas (Richards et al. 2020). The experimental design adopted was the same used in Montagnana et al. (2021): ten landscapes of 5 km radius from the centroids (hereafter regional landscape) with three landscapes of 1 km radius (hereafter local landscape) nested within each regional landscape, totalizing 29 local landscapes sampled. The landscape selection procedure involved searching for an orthogonal gradient of landscape heterogeneity and forest cover and areas with authorization from landowners. All regional landscapes were spatially independent (i.e. did not overlap).
Data collection
In each of the 29 local landscapes, 110 trap-nests (TN) of cardboard tubes (6 cm long and 0.6 cm in diameter) and approximately 20 TN of bamboo segments (several different lengths and internal diameters ranging from 0.5 to 2.2 cm) were installed in forest edges (as nesting is highest in forest edges, Rocha-Filho et al. 2017) located in the center of the landscape. TN of cardboard was inserted into wooden boards and TN of bamboo in PVC pipes, both structures were arranged on iron supports fixed to the ground and covered with acrylic tile (to protect against rainfall).
Due to some logistical limitations, such as the approval of landowners to access forest remnants, we were forced to sample the community in two different years. From September 2015 to March 2016, 15 local landscapes were sampled (Period 1), and from September 2016 to March 2017, another 14 local landscapes were sampled (Period 2). The period from September to March was chosen because nesting activities of solitary bees and wasps are intensified during the hot and rainy season (Camillo and Brescovit 1999; Alves-dos-Santos 2003; Santoni et al. 2008; Gazola and Garófalo 2009). Once a month, TN were inspected with the aid of a flashlight and those occupied, closed with soil or plant materials, were removed and replaced by empty TN. In the laboratory, TNs brought from the field were contained in test tubes (cardboard TN) and PVC tubes (bamboo TN), sealed with porous tape, and kept at room temperature for the emergence of individuals. Emerged individual was frozen, pinned, identified, and incorporated into the Collection of Solitary Bees and Wasps of the Department of Biology, University of São Paulo (USP), Ribeirão Preto, São Paulo. TNs that remained closed for a long time were opened to investigate to the mortality cause.
Community and network metrics
Natural enemy community was described using species richness, abundance (number of emerged adults), parasitism rate (number of attacked brood cells divided by the total amount of brood cells), and taxonomic diversity (∆+). Taxonomic diversity is defined as the average taxonomic path length between any two chosen species according to the Linnean classification of the community sampled (Clarke and Warwick 1998). High ∆+ values indicate a taxonomically diverse community composed of distantly related species, while low ∆+ values indicate a taxonomically more homogeneous community composed of related species (Andersson et al. 2013). Natural enemy species were sampled in 29 local landscapes, and in only two landscapes there was no occurrence of parasitism in TN.
We built a quantitative host-natural enemy network for each local landscape, and to describe their structure we used the robustness metric (R). This metric measures the system’s robustness to losing species from one guild (here hosts) quantifying the number of species from the other guild (here natural enemies) that will become extinct. Host species were randomly removed from the network and an extinction curve was generated by plotting the number of remaining natural enemy species against the cumulative number of host species removed; so, R is defined as the area under the extinction curve (Burgos et al. 2007). R = 1 is a system in which most natural enemy species are not extinct with the removal of most host species, and R = 0 is a system that already collapses with the removal of few host species. Two indices were calculated to measure the trophic specialization of host-natural enemy networks: natural enemy’s niche overlap and linkage density. Niche overlap was calculated using the Horn’s index (Horn 1966), indicating how similar the interactions of natural enemy species concerning the host species they parasitize. The index ranges from 0 to 1, where values close to 1 indicate low trophic specialization and values close to 0 indicate high trophic specialization. Linkage density means whether specialist or generalist natural enemies dominate a network (Bersier et al. 2002). Lower values of linkage density show high trophic specialization and higher values of linkage density show low trophic specialization. Of our 27 networks (in two landscapes there was no parasitism of host nests), 7 had less than 3 host and 3 natural enemy species and for that reason were excluded from all network analyses. Community and network metrics were calculated with vegan v. 2.5-7 (Oksanen et al. 2020) and bipartite v. 2.17 (Dormann et al., 2009) packages, respectively, for R software (R Development Core Team 2012).
Landscape metrics
To calculate the landscape metrics was created a land cover map using satellite images with 1 m spatial resolution available in the ArcGIS® version 10.2 software database. We manually digitalized the landscape feature boundaries at a scale of 1:5000 and classified the polygons into seven classes: river and lake, forest formation, agriculture, forest plantation, pasture, wetlands and, urban area (Fig. 1). Landscape metrics were calculated using the Fragstats software version 4.2.1.603 (McGarigal et al. 2012). To calculate forest cover we used the PLAND metric (% of the landscape covered by the focal class). The landscape heterogeneity metric adopted was the Shannon Landscape Diversity Index (SHDI), which increases as the number of different patch types increases and/or the proportional distribution of the area between patch types becomes more equitable (McGarigal et al. 2012). PLAND and SHDI were calculated at the local and regional scales. CONNECT metric was used to calculate the forest functional connectivity, which calculates the percentage of all possible pairs of forest patches less than 1000 m apart. This distance is recognized as the average flight distance for solitary bees (Zurbuchen et al. 2010). CONNECT was calculated at the local level only.
Statistical analysis
Spatial autocorrelation – Like Osorio et al. (2015), we used the Mantel test to explore the spatial structure of natural enemy species composition. A matrix of geographic distances between the trap-nest locations and a matrix of species composition’s Bray-Curtis quantitative dissimilarity index was used. We found no spatial autocorrelaton (r = -0.02, P = 0.647), so we did not include spatial coordinates in our analyses.
We fitted Generalized Additive Models (GAM) with Gaussian probability distributions. GAM models allow a better interpretation of the results when the relationships between the variables show an absence of linearity (Zuur et al. 2009), and when the same model has two explanatory variables together. The explanatory variables were forest cover (%), landscape heterogeneity at the local and regional levels, and forest connectivity (%) at the local level. Previous studies have shown a positive relationship between the natural enemy community structure and the host community (Albrecht et al. 2007; Holzschuh et al. 2010; Pereira-Peixoto et al. 2016). So, to remove the effect of the host community, we used the residuals as the sole response variable rather than the natural enemy species richness, abundance, and taxonomic diversity (See Appendix A: Table 1); for the parasitism rate we did not use the residuals. For the network analysis, we used as response variables: robustness (R), niche overlap, and linkage density. A null model was also included in the list of competing models, representing the absence of an effect (Martensen et al. 2012).
Models were selected based on the Akaike Information Criterion corrected for small sample sizes (AICc, Burnham and Anderson 2002). The best models were those with the lowest AICc values, and those with AICc values different from the best model but with an ∆AICc lower than 2.0 were considered equally plausible. In cases where more than one model was considered equally good, we weighted the models (wi) using the weight of evidence (wi > 0.5) as the threshold. All analyses were carried out in R Studio software version 1.0.143 (RStudio Team 2022) using the mgcv and bbmle packages (Wood et al. 2011; Bolker and R Development Core Team 2022).