The number of wild bees and cavity-nesting wasps is abundant in agricultural areas and they contribute significantly to ecosystem services. Due to their specialization in nesting sites and  food sources, these groups are sensitive to habitat condition changes and they are therefore important indicators for environmental impact assessments. As semi-natural habitats are steadily declining and often understudied, their significance for research is increasingly recognized. During this research, the role of wild bee species and cavity-nesting Hymenopteran taxa as indicators was examined, along the unique combination of high nature value and traditional land use habitats in Central Europe, Transylvania. Transects and trap nests were used to test the diversity and abundance of wild bees and cavity-nesting Hymenopterans in order to identify possible differences between highly protected and less protected areas. The differences in taxonomic groups between the sites and the potential effects of landscape structure on wild bees and cavity-nesting Hymenopterans were also assessed. We detected a high diversity of wild bee species and a significant species replacement from one study year to another. Among the nest-building Hymenopteran taxa, the majority of nests was built by Trypoxylon sp. during both study years, with a stronger dominance in the second year. The different taxonomic groups of wild bees and cavity-nesting Hymenopterans showed differences in their habitat affinities. The majority of the sampled bumblebee species as well as Trypoxylon sp. had an affinity towards the study sites located within the highly protected study area. Altogether, we found different habitat preferences for different Hymenpteran groups (both wild bees and wasps) and conclude that these groups definitely have the potential to serve as indicators for differences in the intensity of land use. cavity-nesting bees, solitary bees, wild bees, landscape structure, spider prey preference.

Study sites

Our research was carried out in Transylvania (Romania) between 2018 and 2019 in areas with traditional cultivation, extensive farming and high nature value. These areas rich in trees and shrubs have a high plant diversity and density. All study areas were located on the border of the counties of Hargita and Kovászna in a hilly-mountainous landscape, where the average altitude is 580m. These areas are mainly used as mowing fields. Mosaic plots rarely exceed one hectare in size and are relatively far from villages. The mosaic plots are mowed at different times of the year, providing a continuous source of food for pollinators throughout the year. The majority of the study sites was located in the Vargyas valley, which is part of the Natura 2000 network. The importance of these high nature value areas has recently increased due to the intensification of agriculture. There were three study areas, which were located within or in close proximity to Natura 2000 sites (Figures 1a and 1b): The Kormos-valley ('K'; Natura 2000 site ROSCI0091), the northern Vargyas-valley ('NV'; Natura 2000 site ROSPA0027), and the southern Vargyas-valley ('SV'; Natura 2000 site ROSCI0036). The Kormos-valley was the study area closest to human settlements, with small forest patches and a low number of arable fields in addition to extensive meadows. The study area in the northern Vargyas-valley was located at a medium distance from human settlements compared to the other two study areas and is mainly used as mowed pastures, grazed in early spring and late autumn. The nature reserve in the southern Vargyas-valley was located furthest from human settlements and is characterized by forests and mowed areas rich in trees and shrubs.

Sampling of wild bees along a transect

The study was conducted over two years (2018 and 2019), with sampling carried out once in May, twice in June, and once in July using the transect method. In each study area, three sampling sites were randomly assigned. Three pairs of transects of ca. 200m in length were assigned at the sampling sites, where two persons sampled simultaneously (Figure 1c). A distance of ca. 50 m was left between the two transects and the transects did not cross. In all cases, sampling took place between 9 and 17 h, in sunny and calm weather, above 13°C and below 32°C. There was a considerable difference in the temperatures measured during May between the two years (3.7° C; Figure 2a), while the precipitation was very similar (Figure 2b). During sampling, field-identifiable individuals (e.g. honeybee, some bumblebee species) were recorded, and those that could not be identified with certainty were placed in 70% alcohol. The wild bees collected were stored in alcohol and sent to taxonomic expert, Zsolt Józan for identification at the end of each year.

Sampling of cavity-nesting Hymenopterans with trap nests

To record the cavity-nesting Hymenopteran taxa in the three study areas (K, NV and SV), we set up eight study sites in 2018 (the study site NV2 was not sampled in this year) and nine study sites in 2019. Four trap nests in 2018 and three trap nests in 2019 per study site were deployed in 2018 and 2019 (see also Lajos et al., 2021 and Lajos 2022). The trap nests were PVC pipes of ca. 12 cm in diameter and 23 cm in length filled with common reed (Phragmites australis), which were installed in May in both years. The reed was cut with a circular saw to a length of 22 cm, ensuring that each reed contained a node. The trap nests were placed 1.5-2 m above ground on nearby shrubs. The trap nests were collected in early September in both years and then stored in a cold shaded outdoor place until January - February. As the weather warmed, the trap nests were taken from this outdoor place and put into a refrigerator at 4-7°C. At the same time, we started processing the trap nests and recording data. We counted the number of reed stems in the trap nests and then cut each reed stem taking care of the nests inside. Where nests were found they were marked with a unique identification code. The following parameters were recorded for each reed: 1) diameter of the reed, 2) number of nest cells, 3) type of nesting material, 4) colour of larvae or cocoons, 5) type of food found in the brood, whether nectar-pollen mixture or dead arthropods (usually spiders). Seven different nest types were identified based on the parameters 3, 4 and 5. In both years, several reeds per nest type were allowed to grow (minimum 2 reeds in 2018, 15-20 reeds in 2019), which were grown at room temperature. Once the individuals were mature, they were placed in 70% ethanol and then sent for identification. The taxonomist Zsolt Józan identified each two genera of Eumenidae (Ancistrocerus and Symmorphus), Megachilidae (Megachile and Osmia) and Pomilidae (Auplopus and Dipogon) as well as one genus of Colletidae (Hylaeus) and Crabronidae (Trypoxylon). Each nest type was assigned to a genus, except for the two genera Ancistrocerus and Symmorphus of the subfamily Eumeninae, which could not be distinguished from each other based on the nest parameters.

Landscape structure

The landscape structure surrounding the study sites was quantified by calculating the proportion and edge density, which measures the edge length per area of a landscape element and thus measures the fragmentation of these elements, of woodland patches, the majority of which represents the original natural vegetation in this area. Determining the proportion and edge density of a landscape element is in most cases sufficient to quantify the structure and potential effect of this landscape element on different biological variables or processes (Lajos et al. 2020; Lajos et al. 2021a; Lajos et al. 2021b). In order to calculate these two landscape metrics, vector maps with a radius of 250 m were created in QGIS 2.18.9 (QGIS Development Team 2009), around the study sites, using the European Terrestrial Reference System 1989 Lambert Azimuthal Equal Area (ETRS89/ETRS-LAEA; EPSG: 3035) coordinate reference system (Figures 1a and 1b). The centers of these circular vector maps were at an approximately equal distance from the sampling transects for the wild bees and from the positions of the trap nests for the cavity-nesting Hymenopterans (Figure 1c). The vector maps were then transformed into raster maps with a raster cell size of 1 × 1 m, which were used to calculate the landscape metrics ‘Percentage of Landscape’ and ‘Edge Density’ in FRAGSTATS v4.2.1 (McGarigal et al., 2002), applying an eight-cell neighborhood rule.

Data analyses

All statistical analyses were carried out in R v3.6.3 (R Core Team (2020). All graphs were created using the R package ‘ggplot2’ (Wickham, 2016). The species diversity of bees was assessed by calculating the Shannon's Diversity Index (SDI) for each site using the R package ‘vegan’ version 2.5-6. (Oksanen et al., 2019). Non-parametric Wilcoxon signed-rank tests were used to test for significant differences between the abundances and the SDI of the eight most abundant wild bee taxa, as well as between the abundances of nests built by seven different cavity-nesting Hymenopteran per site and study year.

Potential habitat affinities (= towards the study areas and sites, respectively) were investigated with principal component analyses (PCAs) and dendrograms using several different functions from the R packages ‘FactoMineR’ and ‘factoextra’ (Husson et al., 2014). These calculations were based on the numbers of specimens sampled for the eight most abundant wild bee taxa and the numbers of nests built by the seven different cavity-nesting Hymenopteran taxa.

The relationships between the proportion and edge density of woodland patches around the study sites and the numbers of the eight most abundant wild bee taxa or the (occupied) brood cells of the seven different Hymenopteran nest types were analyzed using generalized linear mixed models (GLMMs) assuming a Poisson distribution from the R package ‘lme4’ (Bates et al., 2015), with the study site and year as random effects. The two landscape metrics of woodland were scaled prior to these GLMMs. The residuals of these GLMMs were tested for uniformity, dispersion and outliers using functions from the R package ‘DHARMa’ (Hartig, 2020). We did not find any significant deviations for the residuals of the tested GLMMs in all cases, where significant relationships were detected.

Finally, we also checked for spatial autocorrelation (Moran's I) using the R package ‘ape’ (Paradis and Schliep, 2019) and reanalyzed the data in all cases, where we found significant relationships. Here, we used Generalized Least Squares (GLS) models by REML from the R package ‘nlme’ (Pinheiro et al., 2013), incorporating a Gaussian correlation structure to correct for spatial autocorrelation. The coordinate reference system used for this analysis was ETRS89/ETRS-LAEA (EPSG: 3,035), the same one as used for mapping.