Estimating occupancy and abundance of bees in the mid-Atlantic United States.

Figure 1. Map of study area showing sampling locations (n = 127) in mid-Atlantic U.S..

Daily weather data
Daily precipitation (cm), temperature ^oC, wind speed (m/sec), and humidity (%) were compiled for each site from the closest weather station using the ‘geonames’ (Rowlingson and Rowlingson 2015) and ‘weatherData’ (Narasimhan 2017) packages in R, and then averaged over each sampling week to be included in modeling the detection process in both occupancy and abundance models.

Land cover data
We used 2011 National Land Cover Data [Homer et al. (2015); see Fig. 2], and clipped polygons to study area boundaries. These data comprised of 15 land cover classes, were used to extract relative proportions of each land cover class around each of the sampling locations at 3 spatial scales: 200 m, 500 m, and 1000 m buffer lengths (radius) in R (Team 2016). Land cover data was also used to create predictive habitat-based species occupancy and abundance maps for the mid-Atlantic region.

Elevation data We extracted elevation (m) for each of the sampling locations to be used in occupancy and abundance models.

Figure 2. U.S. Geological Survey 2011 National Land Cover Data used to extract land cover classes.

Figure 3. U.S. Geologial Survey 2011 National Land Cover Data clipped to study area counties featuring 15 land cover classes.

Estimating probability of species occupancy

We analyzed bee count data using R, the open-source statistical computing language (Team 2016). Using the R package ‘unmarked’, we estimated the probability of occupancy for each species using the ‘occu’ function (Fiske et al. 2011). This model is designed for single-season, single-species occupancy estimates that assume population closure during the sampling period. The hierarchical design of the model allows us to account for variables that potentially affect detection probability using a Bernoulli distribution, and is shown in the nested design: \(Zi \sim Bernoulli (ψi)\) where \(Zi\) is the unobserved occupancy state (i.e., present or absent) at site \(i\), and \(ψi\) is the probability of occupancy at site i. This function of a Bernoulli distribution that describes site occupancy driven by species presence or absence at a given location is nested within the model for the detection process: \(yij | Zi \sim Bernoulli (Zi • pij)\) where \(yij\) is the estimated occupancy at site \(i\) and survey \(j\), \(Zi\) is the unobserved true occupancy state (from the function above), and \(pij\) is the detection probability at site \(i\) and survey \(j\) (MacKenzie et al. 2002, Royle and Dorazio 2008).

Estimating species abundance

We estimated species abundance using hierarchical N-mixture models in the package ‘unmarked’ with the function ‘pcount’ (Fiske et al. 2011). Abundance estimates based on raw counts of individuals that do not account for detection probabilities of the individuals counted, can lead to biased population estimates (Royle 2004), and can potentially underestimate population sizes (Royle et al. 2007, Kéry 2008). Hence, we used data from 12 repeated weekly visits to sampling locations to estimate species abundance while accounting for detection probability, and the influence of observation-level covariates (e.g., week, mean precipitation, mean temperature) on detection (Fiske et al. 2011). We used N-mixture models described in Royle (2004) and Kéry et al. (2005) to estimate abundance, which are hierarchical in how the function describing true unobserved species abundance within a given area (e.g., county) is nested within, and in turn influences, the model describing the detection process (i.e., the raw observations). This model assumes closure of the sampled population during the survey period. The model for true unobserved abundance is: \(Nit \sim Poisson (λit)\) Where \(Nit\) is the unobserved true abundance at site \(i\) at time \(t\), which is a function of the Poisson distribution, in which the mean is equal to the variance, with a mean of \(λ\) at site \(i\) at time \(t\). The above model is then nested within and influences the model describing the detection process: \(yijt \sim Binomial (Nit, pijt)\) where \(yijt\) is the observed abundance at site \(i\) during visit \(j\) at time \(t\), which is a function of a binomial random variable with the parameters \(Nit\) (true abundance at at site \(i\) at time \(t\)) and \(pijt\) (detection probability at site \(i\) during visit \(j\) at time \(t\)) (Royle 2004, Kéry et al. 2005).

AIC model selection

We used Akaike’s information criterion (Akaike 1981) to compare \(a priori\) determined models of occupancy and abundance for each species. Models were developed to test for relationships between species probability of site occupancy and abundance and land cover class types: forest, grassland, and developed at varying spatial scales (i.e., 200 m, 500 m, and 1000 m) around each sampling location. Top performing models were then used to estimate predicted values and plotted.

Table 1. List of unique species detected and associated summaries of detections

Figure 4. Site-level averages of raw counts (mean \(\pm\) SE) for 7 selected species. Sites are sorted from higher to lower counts.

Figure 5. Weekly averages for 7 selected species (mean \(\pm\) SE) of raw count data by year. Shown with non-linear model fitted predicted values for each year.

Figure 6. Chloropleth map showing estimated county-level abundance for 7 selected species.

Table 2. Example AIC table for 7 selected species showing model comparisons.

[1] “Andrena erigeniae”

Figure 7. Plot showing relationship between probability of occurence and elevation for 7 selected species.

Figure 8. Predictive map showing probability of occurrence for 7 selected species as a function of elevation (m).

Call: occu(formula = ~week + temp ~ elevation_m, data = umf)

Occupancy: Estimate SE z P(>|z|) (Intercept) 1.23597 0.174731 7.07 1.51e-12 elevation_m -0.00199 0.000962 -2.07 3.89e-02

Detection: Estimate SE z P(>|z|) (Intercept) -0.457 0.1594 -2.87 4.13e-03 week -0.402 0.0408 -9.85 7.04e-23 temp 0.116 0.0251 4.62 3.80e-06

AIC: 2711.771