Decline of six native mason bee species following the arrival of an exotic congener

A potential driver of pollinator declines that has been hypothesized but seldom documented is the introduction of exotic pollinator species. International trade often involves movement of many insect pollinators, especially bees, beyond their natural range. For agricultural purposes or by inadvertent cargo shipment, bee species successfully establishing in new ranges could compete with native bees for food and nesting resources. In the Mid-Atlantic United States, two Asian species of mason bee (Osmia taurus and O. cornifrons) have become recently established. Using pan-trap records from the Mid-Atlantic US, we examined catch abundance of two exotic and six native Osmia species over the span of fifteen years (2003–2017) to estimate abundance changes. All native species showed substantial annual declines, resulting in cumulative catch losses ranging 76–91% since 2003. Exotic species fared much better, with O. cornifrons stable and O. taurus increasing by 800% since 2003. We characterize the areas of niche overlap that may lead to competition between native and exotic species of Osmia, and we discuss how disease spillover and enemy release in this system may result in the patterns we document.


Analysis of intra-season sampling effort variation across time
Variation in sampling effort across years was accounted for in the full mixed models by setting "bowl days" as an offset variable for every cluster-year. However, if variation in sampling effort changed over these years within March to July, it could bias our estimations of increase or decline of catch frequencies based on the different phenology of each species. To examine variation of sampling effort over this scale, sampling events were coded into four time bins covering March to July (the springtime and early summer period in which all Osmia specimens were recorded), with each time bin split as evenly into days as possible from March 1 to July 24, the latest day of any sampling event in the time series (all time bins contained 36 days with the exception of the fourth time bin that contained 38 days). The response variable was the log-transformed sum of all bowl days per year within each time bin. Year, time bin, and their interaction are fixed effects.

Emergence timing among species
We compared relative emergence pattern among the eight most common species in the dataset. Of the 298 spatial-clusters used for overall analyses, there were 147 instances in which a spatial cluster had been sampled at least twice within the same year. Although the haphazard sampling pattern within spatial clusters does not yield certainty regarding the emergence time of any individual species on a given year, when clusters were sampled at least twice in the same year they present an opportunity to see if two species caught in the same cluster were first caught at different time periods. We looked at all species pairs independently, limiting the dataset to all cases in which those two species had been caught in the same cluster and that cluster had been sampled multiple times. We used a Wilcoxon signed rank test to determine if there was a significant difference in first detection between the two species. We used a Bonferroni correction to account for pairwise testing among all species, and refer to significance only for species pairs with a p-value < 0.0017 (Supplementary Table S3).

Landcover types sampled over time
Landcover type of sampling sites was assessed by analyzing percent cover within 200m of sampling sites using the remotely sensed National Land Cover Database 20111 in ArcMap 10.6.12. Landcover types were grouped into five categories: Forest (deciduous, conifer, mixed forest, and woody wetlands); Developed (developed open space, developed low intensity, developed medium intensity, developed high intensity); Open Habitat (shrub/scrub and grassland/herbaceous); Agriculture (pasture/hay and cultivated crops); and Other (all other landcover types). Five simple regressions, one for each landcover type, were performed with percent landcover type as the response variable and "year" as the independent variable.   Table S2. ANOVA output of sampling effort over time. There was a significant difference among years in log-transformed bowl days (as expected and therefore accounted for in the full mixed models as the offset variable), the sampling effort was not significantly different across time bins, and there was no significant interaction effect between year and time bin. The finding of no significant difference in sampling effort across time bins, paired with the lack of a significant interaction between time bin and year gives us reasonable confidence that our analyses are not significantly biased by researcher-based shifts in sampling periods within the months spanning flight seasons of Mid-Atlantic species of Osmia. The failure to find a significant interaction allayed our concern that there was any systematic shift in sampling effort from 2003 through 2017.  Figure S1. Log-transformed sums of bowl days over year, with lines of best fit for each time bin.  Figure S2. Estimated rates of yearly change for models of each species-specific and genus-wide analysis performed, with 95% confidence intervals. Estimates in original scale (reverse i-link) are displayed. Species labeled with an asterisk (*) are exotic to North America. A) results using 3 km clustering; B) results using 5 km clustering; C) results using 10 km clustering; D) results using 20 km clustering. Table S4. Measure of Moran's I of the standardized residuals for each species-specific and genus-wide models used in main analysis. Significance is denoted as **** = < 0.0001; *** = < 0.001; ** = < 0.01; * = < 0.05.

Species
Moran's I  Rate of change (Estimated)

O . b u c e p h a la O . a t r iv e n t r is O . li g n a r ia O . p u m il a O . g e o r g ic a * O . c o r n if r o n s * O . t a u r u s O s m ia (s p p . n > 5 0 ) A ll O s m ia O . c o ll in s ia e
C.  Figure S3. Only multi-sampled clusters used for analyses: xx number of locations. Estimated rates of yearly change for models of each species-specific and genus-wide analysis performed, with 95% confidence intervals. Estimates in original scale (reverse i-link) are displayed. Species labeled with an asterisk (*) are exotic to North America.