Novel approaches to sampling pollinators in whole landscapes: a lesson for landscape-wide biodiversity monitoring

Biodiversity monitoring programs require fast, reliable and cost-effective methods for biodiversity assessment in landscapes. Sampling pollinators across entire landscapes is challenging, as trapping needs to cover many habitat types. We developed and tested a landscape-wide sampling design for pollinators. We assessed the predictability and stability of pollinator biodiversity estimates in agricultural landscapes, and tested how estimates were affected by sampled habitat, landscape composition and spatial scale. We sampled pollinators using pan traps at 250 locations in 10 replicated landscapes measuring 1 × 1 km and calculated bee richness predictions based on different sample sizes. Traps were placed regularly in each landscape, sampling each habitat proportionally to its area. Landscapes contained semi-natural habitats, crop fields and forests and differed in the amount of a mass-flowering crop (oilseed rape). Regular sampling reflected local habitat amount. Compared with cereal fields, significantly more pollinators occurred in oilseed rape, and fewer in forests. Sampling in only one habitat type led to biased estimates of landscape-wide bee species richness, even when sample size was increased. The spatial scale of best predictions depended on the sampled habitat. Species richness was overestimated when sampling was limited to semi-natural habitats and underestimated in oilseed rape fields. Precision increased with the number of sampling points per landscape. To study landscape-wide pollinator biodiversity, we suggest to sample multiple sites per landscape in a broad range of resource-providing habitat types, with sample sizes proportional to habitat amount. Our approach will also be useful for biodiversity monitoring programs in general.


Introduction
In the context of recent reports on declines in pollinator biodiversity (Biesmeijer et al. 2006;Potts et al. 2010) and overall insect biomass (Hallmann et al. 2017;Vogel 2017), the debate has opened up on what the main drivers of these declines could be. For example, the study on insect declines by Hallmann et al. (2017) has received criticism for focusing ''only'' on protected areas, while little is still known on pollinator abundance or richness on cropland or on larger spatial scales. Increasingly, landscapes outside protected areas are moving into focus (Martin et al. 2012;Willis et al. 2012), yet a clear landscape-wide sampling methodology for pollinators is lacking.
While there is a staggering amount of separate studies for pollinators in individual habitat types (e.g. grassland, cropland, forest), almost nothing is yet known on pollinator abundance or even pollinator richness on a landscape-wide scale. Where, in a 1-km 2 landscape, do which pollinators occur? And what happens when major flowering resources in that landscape occur in different amounts (as is the case for example for mass-flowering crops)?
In order to address these questions, and to properly design potential future monitoring programs, an adequate sampling design is needed that covers a large spatial scale (say, square kilometers). In addition, estimates on the required sampling effort are required. A potentially suitable approach involves sampling a range of habitat types on a large spatial scale, for example by (i) selectively sampling several habitat types (Tylianakis et al. 2005;Holzschuh et al. 2016), (ii) establishing transects with their length adapted to local habitat area at nested spatial scales (Gillespie et al. 2017), or (iii) by establishing sampling grids covering several hectares or square kilometers (hereafter termed ''landscape grid method''; Fig. 1).
When two-dimensional maps of pollinator biodiversity patterns across whole landscapes are desired for a realistic estimation of landscape-wide species abundances, a grid-based sampling approach is particularly useful, yet this approach has surprisingly rarely been used so far (Beduschi et al. 2015). Up to now, examples of grid-based sampling at multiple spatial scales come e.g. from a sampling campaign for soil insects, where Benefer et al. (2016) showed detailed maps of the distribution of various taxa on sites measuring about 1 9 1 km. Similar approaches were used in the pan-European project ''Greenveins'' for insects and birds (Dormann et al. 2007b;Le Féon et al. 2010), and grid-based approaches in general are widely employed in biodiversity monitoring schemes (e.g. Manley et al. 2004) or forest inventories (Kowalski et al. 2011).
Establishing sampling grids (grid-based sampling, also termed ''regular '' or ''centric systematic sampling'';Krebs 1999;Ripley 2005) allows samples to be taken evenly across a landscape. According to Ripley (2005), systematic sampling is ''best unless [there is] a strong periodicity'' in the response variable that coincides with the sampling interval. A major prerequisite for grid-based sampling (of pollinators) is that each major habitat in a landscape should be sampled at least once (termed ''spatial lag''; Fortin and Dale 2005).
In the present study, we use a design where samples are placed regularly throughout the landscape (Fig. 1). By imposing a regular sampling grid on a given landscape, habitats are sampled proportionally to their area in the landscape. Our aim is to investigate how differences in the number of samples taken influence estimates of pollinator species richness. We sampled pollinators in 10 replicated landscapes, each comprising 25 sampling points. These landscapes had been intentionally selected a priori to differ in the amount of a mass-flowering crop (oilseed rape-Brassica napus L., also termed canola). Proportion of oilseed rape (OSR) was chosen as a likely determinant of pollinator species richness, as mass-flowering crops are known to have strong effects on pollinator biodiversity (e.g. Westphal et al. 2003;Diekotter et al. 2010;Holzschuh et al. 2016). Of course, we could have used another gradient (e.g. in amount of arable land or habitat connectivity), but in our case OSR was known to be an important crop attracting large quantities of pollinators. We therefore knew a priori that OSR would be a strong explanatory variable. We tested the following hypotheses: (1) Grid-based sampling allows to estimate pollinator biodiversity on a landscape scale, and pollinator richness for a wide range of habitat types can be determined. This is because a regular sampling grid will correlate with landscape-wide abundance of each habitat. (2) The proportion of oilseed rape at different spatial scales (e.g. 0-100 m, 250-500 m etc.) around each sampling point will affect estimates of pollinator biodiversity. This is because pollinators respond positively to the abundance of mass-flowering crops in a landscape (as long as it is flowering). (3) Sampling in all habitats, only in OSR fields or only in seminatural habitats affects pollinator biodiversity estimates, because of habitat-specific species compositions.

Assessment of pollinator biodiversity
The study was performed in 10 landscapes in the surroundings of Göttingen (51°32 0 N, 9°56 0 E) in Central Germany in 2011 (Fig. 1). The landscapes measured approximately 1 km 9 1 km (mean area ± SD: 0.93 ± 0.23 km 2 ) and represented a gradient of percent area occupied by oilseed rape fields. Care was taken that all other habitat types (grassland, forest, cereal fields, root crops, corn fields) were present in each landscape. Sites were selected a priori out of a total of 13 potential sites that had been visited in the field, ensuring that OSR was statistically independent of amount of other habitat types, and sites with low or high OSR were spatially interspersed. In each landscape, sampling was performed in a grid measuring 1x1 km, comprising 5 9 5 points ( Fig. 1), which was laid out over the landscapes to always include forest margins and grasslands (semi-natural habitats) as well as crop fields, while excluding cities or villages. 1 9 1-km grids were first roughly placed in Google Earth (Ó Google, Inc.), and final position was decided after extensive field visits. Two areas were slightly smaller than 1 km 2 to exclude settlements. c Detail of one of the ten sampling grids. The x and y axes are the coordinates in the World Geodetic System (WGS84) with a transverse mercator projection (latitude of origin = 0, longitude of origin = 9, scale factor = 1). a Grid points in a 1 9 1-km landscape between the villages of Barlissen and Atzenhausen (Southwest of Göttingen, Germany). Points are medoids of three independent handheld GPS measurements. b interpolated map showing pollinator species richness (high richness indicated by lighter colours). Interpolation was done using gridded bivariate spline interpolation for irregular. (Color figure online) Sampling habitats included oilseed rape fields, cereal fields and semi-natural habitats, which comprised grasslands and forest margins. Satellite-based image classification was used to assess landscape composition (OSR, cereals, grassland, forest, corn, root crops; see Fig. 1b) for each landscape (full 1 9 1km scale), and at six additional spatial scales. These scales were represented by six nested rings (sensu Schneider et al. 2011) with the following radii: 0-100 m, 100-250 m, 250-500 m, 500-750 m, 750-1000 m and 1000-1500 m, using ESRI Ò Arc-Map TM 10. False-colour satellite imagery was provided by RapidEye TM . Image classification was performed using ENVI EX (ITT Visual Information Solutions GmbH, Gilching, Germany), using a training dataset (ATKIS DLM 25/1, ATKIS 2010).
Yellow pan traps filled with water (0.75 L volume, 156 mm diameter) were mounted to a wooden pole approximately in the center of each cell of the grid (centric systematic sampling method) and exposed for 3 days in May and June 2011; however, we focus here on the June dataset (after oilseed rape flowering), because the highest bee species richness can be observed during this month (Holzschuh et al. 2011). Traps were placed at vegetation height and exactly at the location given by the grid-based sampling scheme, but avoiding roads or farm tracks. As three pan traps were damaged, we had 247 samples overall. All wild bees were sent to specialists for identification.

Data resampling
To determine how the number of samples per landscape affected the results, we randomly sampled 5, 10, 15 and 20 points per landscape from the full dataset (N = 247), with each of the new datasets subsequently analyzed (N = 50, 100, 150 and 200, respectively). This was repeated 50 times for each subset of number of points, resulting in four sets of model results each with 50 outcomes.
To assess the effects of differences in the locally sampled habitat for each landscape, we took two subsets of our data that included samples collected only in semi-natural habitats (five points per landscape; N = 50) or only in oilseed rape fields (mean points per landscape ± SD 6.2 ± 5.65; N = 56). These habitats were chosen as they represented two extremes in land-use intensity: The semi-natural habitats are often of high conservation value and tend to be preferentially sampled in ecological studies. The oilseed rape fields represent very homogeneous agricultural areas. In each landscape, one point within the chosen habitat was sampled at random, creating a new dataset (N = 10) that was subsequently analyzed (see section Statistical Analyses). This procedure was repeated 50 times per subset (semi-natural habitats and oilseed rape fields) to obtain a wide range of possible results, yielding two sets of model outputs, one for semi-natural habitats and one for oilseed rape fields, each containing 50 outcomes.
Overall, this resulted in the following three datasets used for statistical analysis: (1) all data points collected following the landscape-grid approach.
(2) all data points collected only in semi-natural habitats (habitat-selection method); (3) all data points collected only in oilseed rape fields (habitat-selection method); and A summary of the resampling methods can be seen in Table 1. The selection of points for the new datasets was always repeated 50 times, with each of these new datasets analyzed accordingly. The complete datasets were also analyzed to detect the effect of sampling only one kind of habitat several times per landscape.

Statistical analyses
All analyses were performed using R 3.5.1 (R Core Team 2018). Overall effects of local habitat on pollinator species richness were analyzed using mixed-effects models with ''landscape'' as a random effect, and an exponential variance function to account for heteroscedasticity. Local habitat was the only fixed-effects term in these models.
We determined the relevant spatial scale(s) using linear models fit by generalized least squares (GLS) as these models allow an explicit incorporation of spatial autocorrelation by fitting a variance-covariance matrix (Dormann et al. 2007a). As generalized least squares models returned a log-likelihood, this allowed us to use information theoretic approaches (AICc) for model selection. The response variable was bee species richness and was log transformed (ln (y ? 1)) to restrict predicted values to be nonnegative. The initial explanatory variables were the proportions of the area occupied by oilseed rape within the six spatial scales. All models were simplified using a modified version of the stepAIC function (Venables and Ripley 2002), corrected for sample size (AICc; Burnham and Anderson 2004); note that the penalty term in AICc converges to zero for large N. When more than one point per landscape was sampled, we defined a spherical correlation structure using the coordinates of the sampling points and the landscapes as a grouping variable to account for spatial autocorrelation. GLS models were fitted using the function gls from the ''nlme'' package 3.1-137 (Pinheiro et al. 2018). For maps showing actual sampling locations (Fig. 1a-c), we calculated averages of three measured GPS coordinates by clustering of input data around N = 25 medoids per landscape (R function ''pam'' in package ''cluster''; Maechler et al. 2017).

Results
Using grid-based sampling allowed us to sample habitats according to their true amounts in the landscapes (Fig. 2): The total area covered (Fig. 2a) was almost perfectly matched by the number of sampling points (''selection frequency'' in Fig. 2b), indicating that our sampling methodology indeed reflected landscape-wide habitat amounts.
Overall, we collected 76 bee species, excluding Apis mellifera (Linnaeus, 1758). Thirty per cent of the species (23 spp.) were not found in semi-natural habitats and 55% (42 spp.) were not found in oilseed rape fields. Across all landscapes, bee species richness was significantly higher in OSR fields than in cereal fields, and significantly lower in forest habitats than in cereal fields (Fig. 3). Interpolation allowed landscapewide prediction of pollinator species richness (Fig. 1c).
The number of sampling points had large effects on the detection of relevant spatial scales in models on pollinator richness vs. OSR among (Fig. 4): When only one point per landscape was sampled (total N = 10), often no scale (i.e. radius of landscape sector) was selected as relevant (72% of the times for seminatural habitats and 56% for oilseed rape fields; Fig. 4a) or found to be significant (80% and 60% of the times for semi-natural habitats oilseed rape fields, respectively; Fig. 4b). Additionally, no clear pattern  identifying a preferred radius was recognized. However, with an increasing number of points sampled in a landscape, it was possible to detect an increasing precision, with the 750-1000 m scale chosen as relevant and significant in the majority of the models (Fig. 4a, b). When 20 points per landscape were sampled (total N = 200), this scale was statistically significant in 94% of the 50 models performed (Fig. 4b).
When all samples collected in oilseed rape fields (N = 56) were analyzed in a single model, only the 1000-1500 m radius remained in minimal adequate models (Table 2). This means that this is the only scale that can explain the data. The model using the dataset restricted to semi-natural habitats (N = 50) selected the same scale as the model that included all N = 247 systematically sampled points (750-1000 m). Nonetheless, the estimate from this full semi-natural habitat model was very different from the one incorporating all habitats (Table 2; Fig. 5). In the semi-natural habitat model, none of the points sampled presented a proportion of oilseed rape greater than 0.3 and only a few exceeded 0.2. This constitutes a truncated oilseed rape gradient, which means that part of the range of the environmental variable was not included in the sampling frame (Albert et al. 2010). As a result, the expected number of bee species in the missing range was clearly underestimated in the outcome of the model, when compared to sampling all habitats.
The estimates of bee richness in relation to percentage of oilseed rape fields, when considering only one point per landscape (N = 10), were very variable, independent of the sampling habitat selected, and fluctuated from negative to positive values (Fig. 6). Furthermore, we found a gradual increase in the precision of the estimates and reduction in bias with a growing number of points included in the sampling, as a larger proportion of the models approached the estimate of the complete model including all N = 247 points sampled (Fig. 6).

Discussion
Our study shows that pollinator species richness can be assessed on a landscape-wide scale using regular grid-based sampling. The sampled habitats closely reflected the ''true'' amounts of habitats in the landscape. The number of samples per area and sampling habitat affect the estimation of landscapewide pollinator species richness.
We found that limiting sampling to only one point (one habitat patch) per landscape yields biased estimates, given that individual points are subject to local stochasticity. Estimates depended on the sampling points chosen, as a consequence of the great variation found among possible sampling points in the landscape. If only a few points in a landscape are sampled, sampling intensity may be increased by exposing traps for a longer period of time (Beduschi et al. 2018).
Under low sampling intensity per landscape, all of the considered radii had equal chances of explaining the data (Fig. 3); thus, the relevant spatial scales in a dataset can only be reliably estimated if either (i) the precision of individual estimates per landscape (i.e. the number of sampling points) is increased or (ii) the number of landscapes sampled is increased (which was not possible in our case). Our simulation study showed that the probability to identify the most influential (''correct'') radius for a given response variable increases with sample size. This sheds new light on previous studies that extrapolated to the landscape scale, but used only a few points in the landscape (not arranged in a grid): For example, several studies have focused on the effect of proportion of oilseed rape fields in the landscape on pollen beetles Meligethes aeneus (Fabricius, 1775), a pest of oilseed rape, reaching very different conclusions. Our results show that limiting the sampling to one habitat type can lead to biased estimates in that they cannot be extrapolated to the whole landscape. This was observed even when the number of samples in that habitat was increased. This can happen, as was the case with the semi-natural habitat samples, because ecological studies often do not encompass the full range of possible environmental conditions (so-called ''truncated gradient''). Truncated gradients can be avoided by sampling across a wider range of environmental conditions (Mohler 1983), as is the case in grid-based sampling.
The implications of our findings are potentially farreaching, as they may also affect species distribution models; if these models are parameterized based on   Edwards et al. (2006) compared how predictions based on a probabilistic versus a non-probabilistic sampling designs reflect the real pattern of lichen species distribution, and found that a systematic grid sampling produces more realistic results than a purposive sampling strategy (where sampling effort is locally increased to detect particular species).
One of the most prominent examples of purposive sampling of pollinators (and other flying insects) is a recent study on insect declines in protected areas (Hallmann et al. 2017). The study reported that insect declines were ''independent of land use composition at surroundings'' of their study locations, but sampling was restricted to protected areas of low productivity. Thus, conservation decisions on a landscape scale, derived from purposive sampling, have to be handled with care, and more studies encompassing a wider range of habitats in a landscape are needed.
It is not surprising that sample size was important for the precision of our estimates. This has already been pointed out by Hirzel and Guisan (2002), while Albert et al. (2010) argued that sampling design and not sample size is the most relevant factor influencing parameter estimation. Additionally, Marsh and Ewers (2013) claimed that both configuration and number of sampling points affect beta-diversity estimates, which results in incorrect diversity partitioning estimates. And even though the present study focused on alpha diversity, we observed that both sample size and sampling design play a significant role, influencing precision and bias (for beta diversity, see Beduschi et al. 2018). Therefore, it is advisable to sample study areas multiple times to reduce uncertainty around the estimates. While this procedure can generate spatial autocorrelation in the residuals, a variety of statistical methods can be used to account for this (e.g. Dormann et al. 2007a). Relationships between bee species richness (log transformed) and proportion of area occupied by oilseed rape within a buffer area ranging from 750 to 1000 m distance from the sampling point. Yellow (dotted), green (dashed) and black lines show the predictions made by generalized least squares models for all the data points collected in oilseed rape fields, seminatural habitats and following a grid throughout the landscape, respectively, as seen in Fig. 5. Each grey line represents the outcome of a generalized least squares model performed in each of 50 datasets created according to the following rules: a 1 point per area sampled in semi-natural habitats; b 1 point per area sampled in oilseed rape fields; c 5 random points per landscape (top right); d 10 random points per landscape; e 15 random points per landscape (bottom centre) and f 20 random points per landscape Finally, we found that the spatial scale determining pollinator species richness also changes with sampling habitat. This indicates that the processes affecting pollinator diversity operate at different scales according to habitat type. For example, the radius of a landscape sector best predicting species richness in oilseed rape fields was larger than in semi-natural habitats, which indicates that landscape-scale dilution effects take place at larger scales as bees spillover to farther areas. This shows that studies that sample only one habitat are valuable to determine how diversity relates to environmental variables or how it increases with area within that habitat type. Nonetheless, it should remain clear that the results will only allow indirect estimates of how the surrounding landscape affects local diversity.
It should be noted that the aim of our study was not to justify pan-trap sampling over other sampling methods. However, we were surprised to see that such a comparatively simple method can produce pollinator richness estimates across a wide range of habitats (as has also been shown, e.g., by Westphal et al. 2008). Future studies could employ colorless pan traps (e.g. Everwand et al. 2014) to avoid too attractive traps.

Conclusions
We demonstrated that the sampling design can affect the predictability of landscape-wide pollinator biodiversity estimates. Our results show that number of samples per study area affected the precision of parameter estimation and the preferential selection of habitats for sampling generated biased estimates of parameter and species richness. Parameter estimates obtained by sampling in only one habitat type may be relevant when the researcher aims to understand biological responses within the boundaries of the habitat. However, studies performed in only a single habitat type cannot be extrapolated to the whole landscape, which is the scale driving population dynamics including extinction and survival, and should therefore be interpreted cautiously. For studies attempting to understand how pollinators respond to landscape components, we suggest that the range of the sampling area, variety of sampling habitats and the number of sampling units should be increased to all habitat types of the landscape level to obtain more reliable results.
We showed how differences in sampling design affect responses of pollinators to landscape-level variables. However, our findings have wider implications also for other organisms with different foraging ranges or home range sizes. Spatially interpolated maps of organism presence should be based on random or regular sampling designs, and grid-based sampling schemes should be employed for biodiversity monitoring schemes as a whole. Some countries, such as Switzerland or Austria, have already implemented grid-based sampling schemes, and it is hoped that the present manuscript will contribute to improving such sampling regimes across taxa.