Hyperabundant black-tailed deer impact endangered Garry oak ecosystem floral and

Garry oak ( Quercus garryana ) and associated ecosystems (GOAE) are culturally and ecologically significant landscapes home to over 100 Species at Risk. With less than 5% of their original extent remaining in Canada, GOAE require urgent attention to guide management of remnant habitats and protect them from ongoing cumulative threats. The loss of top predators has played a large role in a trophic cascade where hyperabundant black-tailed deer ( Odocoileus hemionus columbianus ) are reshaping plant and animal communities, even within protected areas. The impacts of deer browsing on plant communities is well-documented; however, the impact on pollinators that are essential for the restoration and conservation of GOAEs is still unknown. We tested the hypothesis that deer browsing has an indirect impact on bumblebee populations through their direct negative impact on floral communities. We surveyed deer density, floral communities, and bumblebee abundances across ten islands with differing levels of deer presence in the Salish Sea. We found that deer had a negative effect on the abundance of flowering plants, with an even greater effect on native species that bumblebees rely on for foraging. Deer presence had an indirect negative impact on female bumblebee abundances through the depletion of floral resources. Deer also had a negative impact on male bumblebee abundances, potentially caused by lowered colony success over the season on islands with less floral resources. These findings highlight the importance of immediate deer management to prevent the loss of native pollinators and further, potentially irreversible, ecological degradation of GOAEs.


Introduction
Globally, anthropogenic disruption of community and food web dynamics are leading to the degradation of native ecosystems (Estes et al., 2011;Ripple et al., 2014;Wolf and Ripple, 2017).In particular, human removal of top predators has contributed to trophic cascades where unregulated browsing by hyperabundant herbivores has caused degradation of ecosystem function and resilience, leaving ecosystems more vulnerable to large-scale shifts in community composition (Atkins et al., 2019;Beschta and Ripple, 2009).To restore these ecosystems effectively, an understanding of both historic and current ecosystem dynamics is needed (Dudley et al., 2020;Munawar et al., 2020).
Garry oak (Quercus garryana) and associated ecosystems (GOAE) are some of the most at-risk ecosystems in North America, and have undergone dramatic loss and degradation since European colonization (Barlow et al., 2021;Fuchs, 2001).Once widespread throughout Western North America, GOAEs have decreased by a devastating 95% in Canada over the past two centuries (MacDougall et al., 2004).GOAEs contain a plethora of unique native species, culturally significant food plants (Turner, 2020), and over 100 species listed as endangered, threatened or vulnerable by the Committee on the Status of Endangered Wildlife in Canada (COSEWIC) (MacDougall et al., 2004).Historically, Garry oak ecosystems were stewarded by First Nations through cultural burning, harvest of large herbivores such as deer, and the planting, propagating, and harvesting of certain plant species in what are known as root gardens (Martin et al., 2011;McDadi and Hebda, 2008;Turner et al., 2013).GOAEs are considered ecocultural landscapes as a result of their tightly coupled relationship with First Nations stewardship for millennia (Pellatt and Gedalof, 2014).The arrival of colonizers in the region in the mid 18 th century led to the loss and alienation of First Nation communities from their lands putting an end to many of these land stewardship practices.At the same time, top predatorsincluding wolves (Canis lupis), cougars (Puma concolor), and bears (Ursus americanis)were viewed as a threat to settlers and their livestock and extirpated (Arnett, 2000).With the loss of local predators and traditional practices to keep herbivore populations in check, along with a decline in hunting in recent decades, native black-tailed deer populations (Odocoileus hemionus columbianus) have burgeoned up to 10 times their historical abundances with cascading impacts on native vegetation communities (Arcese et al., 2014;Gonzales and Arcese, 2008;Martin et al., 2011).
Remaining GOAEs in Canada are found on south-eastern Vancouver Island and nearby Gulf Islands within the Salish Sea.The hyperabundance of black-tailed deer in the Gulf Island archipelagos of Western British Columbia has been well documented as a major contributor to ecological degradation within the region (Arcese et al., 2014;Martin et al., 2011).Deer browsing has dramatically altered the structure, biomass, and composition of understory vegetation, and has specifically reduced the abundance and presence of ecologically and culturally important native plant species, such as camas (Camassia quamash, C. leichtlinii) and seablush (Plectritis congesta) (Gonzales and Arcese, 2008).Loss of understory vegetation has in turn led to significant negative impacts on understory-dependant songbirds, including the rufous hummingbird (Selasphorus rufus) and Bewick's wren (Thryomanes bewickii) (Martin et al., 2011).We predict that loss of understory floral resources due to hyperabundant deer populations within GOAEs may also be leading to a decline in native pollinators, specifically bumblebees.Although deer hyperabundance is a wide-spread problem across North America, few studies have investigated how wild, predator-free ungulate populations are impacting native pollinator communities (Black et al., 2011;DeBano et al., 2016).
Native pollinators are an important component of most ecosystems and are becoming increasingly recognized for their value in both natural and agricultural settings.Bumblebees (Bombus sp.) are generalist foragers that are integral to the success of native plant communities, and are considered one of the most important taxa of native pollinators (Cameron et al., 2011;Cameron and Sadd, 2020).In the past two decades there has been a dramatic decrease in many bumblebee species in North America, Europe, J o u r n a l P r e -p r o o f and South America (Aizen et al., 2019;Cameron et al., 2011).The abundances of several common bumblebee species have declined by up to 96% in North America, with habitat degradation and loss being identified as one of the major drivers (Cameron et al., 2011).Bumblebees rely on diverse floral cover to obtain necessary proteins, lipids, essential amino acids, and sterols for individual and colony survival and health (Brunner et al., 2014;Vaudo et al., 2016).With degraded habitat and increased fragmentation, obtaining high quality floral resources can become a challenge for bumblebees, leading to increased energy costs in foraging and nutritional stress in the colony.Increased fragmentation also can lead to greater mortality of workers as they have to increase their foraging distances and lengths to obtain sufficient resources, exposing them to more risks.
In this study we investigate the impacts black-tailed deer have on floral and bumblebee communities in Garry oak and associated ecosystems.The Gulf Islands present a perfect study system as they are made up of hundreds of islands and small islets with a wide range of flora and deer density.Although bumblebees play a key role in pollination within GOAEs (Rammell et al., 2018), little is known about their current abundances or distributions within the Gulf Islands.This region is home to the SARA (Species At Risk Act) listed threatened Western bumblebee (B.occidentalis), and this study also assesses current levels of this once wide-spread species (Canadian Wildlife Service and Committee on the Status of Endangered Wildlife in Canada, 2014) on selected islands.Our study examines the floral and bumblebee communities within GOAEs on 10 different islands with varied levels of deer density, and compares floral and bumblebee abundances between sites.Many of these islands are under federal and provincial management, and are under ongoing restoration efforts.Our study provides direct insight into how current management practices are serving floral and bee communities and will help inform future conservation action across the Salish Sea and in other ecosystems facing similar browsing threats across North America (Barlow et al., 2021).

Study sites
This study was conducted in the Gulf Islands of the Salish Sea, off the western coast of British Columbia, Canada (Figure 1).The area lies within the Coastal-Douglas fir Biogeoclimatic zone, which has a mild and temperate climate (Meidinger and Pojar, 1991) with average temperatures ranging from 4 -8ºC during the winter and 12 -18ºC during the summer (Government of Canada, 2021).However, summer 2021 experienced an abnormally hot period ('heat dome') occurring between June 20 th -29 th , with temperatures rising to 40ºC on Vancouver Island, impacting the vitality of many of the plant and animal communities of the islands (Overland, 2021).
Fieldwork was conducted throughout early spring to mid-summer 2021, with 6 survey rounds completed between April and July.There were a total of ten 1-ha study sites, with four sites on islands without deer (Owl, Halibut/SISȻENEM, Mandarte, and Portland) and six sites on islands with deer (Rum, Wallace, D'Arcy, Tumbo, Prevost, and Salt Spring) (Figure 1).These islands were selected for their similar range in size (excluding Salt Spring) and habitat type.When selecting islands, efforts were made to minimize differences in human development and distances to source islands.However, due to restricted access to certain islands, some differences were inevitable.To account for these differences, island size and isolation (from the closest large source island, either Vancouver Island or Salt Spring Island) were included as explanatory variables in model analyses.Although Salt Spring Island differs in size and development from the other islands, it was included because it is typical of developed, human inhabited islands within the Salish Sea and therefor provides an important reference point for discussing deer management on large inhabited islands, not just small islands.Study sites were situated in Garry oak and associated ecosystems which comprise savannas, coastal bluffs, and maritime meadows, and all 1 ha sites included an approximately equal mix of these ecosystem types.Total tree stand basal area was measured J o u r n a l P r e -p r o o f in each study site using the Bitterlich stick sampling method to account for potential differences in tree density among sites (Mulyana et al., 2018;Wilson, 2011).

Deer Density
Deer density was estimated through faecal standing crop (FSC) calculations.FSC was measured by two people slowly walking on either side of a 100m transect, randomly stratified within Garry oak ecosystem habitat of each study site.While walking along the transect both observers counted the number of faecal clumps (defined as having 10 pellets or more) within 1 meter either side of the transect tape, yielding a total search area of 200 m 2 (Martin et al., 2011).Three separate FSC transect walks were done throughout each study site with deer present.These values were input into the following equation to estimate deer density and then averaged (Martin et al., 2011): Where A = Transect walk pellet-clump count (FSC) B = Known deer density of reference island C = Pellet-clump count of 200 m 2 transect walk on reference island This relationship was determined using known deer density values collected from Piers island, where B= 0.22 and C = 5.3 (Martin et al., 2011).
J o u r n a l P r e -p r o o f 2.3 Floral richness and abundance Floral abundances at each site were measured in twenty-five 1 m 2 quadrats.Quadrats were placed across each of the 10 study sites, randomly stratified within Garry oak ecosystem habitat of each study site.In each quadrat, the number of flowering forb and grass stems were counted by species and the percent cover of bare area was approximated (Wray et al., 2014).For large shrubs that had their main stems originate within the quadrats we counted the number of major flowering branches, which were classified as all major branches growing from the main stem.If a plant species could not be identified in the field, a sample specimen (root, stem, and flower) was collected for later identification in the lab.

Bumblebee richness and abundance
The abundance of each bumblebee species was counted by an area search within the 1 ha study sites.One individual slowly circulated through each 1 ha site, moving from flower patch to flower patch for 30 minutes, and all female bumblebees observed visiting a flower (defined as actively touching a reproductive part of a flower) were identified to species (Gillespie and Elle, 2018).Male bumblebees were not identified to species as they are difficult to identify in the field.The species and sex of each bee, the plant species the bee was seen foraging on, and whether the bee was observed foraging on a flowering plant above 1.5 m in height were all noted.Above 1.5 m, deer are generally unable to reach, therefor this threshold is useful for delineating the browsing impact zone.Some of the smaller islands (Owl and Halibut) do not have a contiguous 1 ha area of GOAE, so multiple partial sites were surveyed together to achieve a total survey area of 1 ha, with the timer being stopped when moving from one portion to the other.All bee surveys were conducted by the same individual (KB) to ensure consistency and reduce observer error.If any bumblebee could not be identified on-wing during the survey, it was caught by net for closer inspection.The timer was stopped as netted bees were handled and processed (Boyle et al., 2018;Wray et al., 2014).The majority of caught bees were put in a vial and placed in a cooler until being pinned after returning from the field that day.All netted queens and worker bees that were clearly identified in the vial without the need for further inspection under a microscope were released in the same site they were found.Approximately 15 bees were kept from each survey round as voucher specimens, with an attempt to have equal numbers of all surveyed species.Bumblebee surveys were conducted between the hours of 9:30 and 16:30, during peak bee activity, and limited to sunny days with at least 12ºC and a wind-speed no higher than 3.5 m/sec.Time of surveys were rotated with each site visit, ensuring each site was visited in the morning, mid-day, and afternoon in alternating order throughout the season.These conditions vary slightly from the literature (Gillespie and Elle, 2018); however, the surveyed islands are exposed and experience higher wind speeds compared to inland studies in these areas on average.Hourly temperature and wind levels were recorded during each bumblebee survey from Windy (a mobile app using our current location).

Deer & Flora Relationships
A Linear Mixed Model (LMM) was used to test if deer presence had a relationship with the stem count of all surveyed plant species.A second LMM was used to test the relationship between deer presence and stem counts of only 'usable' plant species, in order to analyse the direct impact on plant species important for pollinator foraging.Plant species were classified as 'usable' by bumblebees through a combination of personal observations of flower visit behaviour, unpublished colleagues' field work, and data compiled in the Pollinators of British Columbia shiny app ( https://shiny.rcg.sfu.ca/u/lmguzman/bc_bees/,Guzman et al., 2020), and erred on the side of inclusion, where if any source showed bumblebees visiting a plant species it was included as 'usable' (Supplemental Materials, Table S1).Both LMMs were run using the lme4 R package (Bates et al., 2015), using deer presence as a fixed effect, and site as a random effect.The J o u r n a l P r e -p r o o f model structure was chosen from an Akaike information criterion (AIC) ranked output of reduced models (Table S2) using the 'dredge' function (MuMln R package) to compare all possible combinations of all fixed-effect covariates, which included deer presence, deer density, average tree basal area, island size, and island isolation (Barton and Anderson, 2020).In the LMM that included only 'usable' flowers, plants with flowers higher than 1.5 meters above the ground were not examined as they represent resources that could not be browsed by deer.For both LMMs, the log of average stem count was used as the counts were on an exponential scale.
Differences between plant species assemblages on islands with and without deer were assessed using Non-Metric Multidimensional Scaling (MDS) ordinations with the Bray-Curtis dissimilarity measure.The position of sites on the ordination relative to other sites is based on the plant species compositional similarity of those sites, where sites with a similar plant composition are close to one another in ordination space, whereas sites with contrasting plant species composition are located further apart.Analysis of Similarities (ANOSIM) was used to test whether the variation in plant species composition between site groups (browsed or unbrowsed) was greater than or equal to the similarity within site groups.Two MDSs and ANOSIMs were run, one comparing total surveyed floral communities and a second comparing only 'usable' floral communities below 1.5 m (accessible to deer browsing).Analyses were completed in R Studio with the Vegan package, using the 'metaMDS', 'vegdist', and 'anosim' functions (Oksanen et al., 2020;R Core Team, 2020).Native plant species richness of usable plants under 1.5 m from the entire season was also calculated for browsed and unbrowsed islands and compared using a 2sided Welch t-test to account for the differences in sample size.

Bumblebee Models
Total bumblebee abundance was analyzed using a Generalized Linear Mixed Model (GLMM) framework, using negative binomial as the family link to account for any overdispersion of the bee count data.Model selection was done in a similar way to the LMMs using the dredge function to determine which variables to include (Supplemental Materials, Table S3).The fixed covariates included in the model were mean site tree basal area, deer density, deer presence (binary), island isolation distance, island area size, mean 'usable' stem count, temperature, time of day, and wind, along with site as the random effect.The top ranked models were tested for uniformity and dispersion using the package DHARMa (Hartig, 2021).Only bees observed foraging below 1.5m were included in these analyses to assess the direct floral use impacted by deer browsing, as well as to reduce identification error that may occur with bees further away that are harder to see or catch.Analyses were completed in R Studio with the lme4 package for the GLMMs and the MASS package for theta approximations used in the GLMMs (Bates et al., 2015, p. 4;R Core Team, 2020;Venables & Ripley, 2002).
For all analyses, male and female bumblebees were analyzed separately, as males and females had different periods of emergence throughout the season.Very few females were observed after July 21 st (the final day of the fifth survey) so only surveys 1 through 5 were used for the female analyses.Likewise, no males were observed before the fourth survey, so surveys 4 through 6 were used for the male analyses (Figure S1).Female analyses include both emerging spring queens and their subsequent workers.Field work was completed before the emergence of new queens in late summer, and thus they were not included in these analyses.MDS was performed to assess the difference between bee assemblages on islands with and without deer and ANOSIM was applied to test if bee species composition varies more between site groups (browsed or unbrowsed) or not (Martin and McIntyre, 2007).

Piece-wise Structural Equation Modelling
To evaluate the direct and indirect relationships of deer presence on floral and bumblebee abundance a piece-wise structural equation model (PSEM) was employed.A PSEM analysis was selected because deer J o u r n a l P r e -p r o o f variables were not included in the female bumblebee GLMM due to AIC dredge ranking (which we interpret to mean deer do not have a direct impact on bumblebees).However, by combining separate models that include deer and bumblebee variables with a PSEM we can determine if deer have an indirect impact on bumblebees instead.Using the piecewiseSEM package, two model equations were included in the analysis: the female bee abundance GLMM, and the LMM of 'usable' stem count, which both met the minimum sample size assumption (Lefcheck, J.S., 2016).The fit of the PSEM to the data was tested using Fisher's C test, as well as a comparison of R 2 values which are more indicative of a good fit in a small two-part PSEM than a Fisher's C test (Shipley, 2000).

Deer density & Floral diversity
Site-specific deer density, as estimated from FSC density within each 1 ha study site, ranged from 0.75 to 1.69 deer per hectare in GOAE on islands with deer present year-round, as calibrated from the reference island with known deer densities (Table 1).
Table 1.Island GOAE sites with their respective deer densities and standard deviations, 'usable' native plant species richness of the entire season, island size and isolation from Vancouver Island or Salt Spring Island.Deer densities were calculated using pellet counts calibrated using known deer density and pellet density on Piers island calibration (see Martin et al., 2011).Deer density values were averaged from n=3 transect walks.The linear mixed models found that deer presence had a significant negative relationship with abundance of all surveyed flowering stems (deer presence: beta= -1.19, SE=0.49,P=0.02), and a slightly greater effect on the abundance of 'usable' flowering stems below 1.5 m (deer presence: beta= -1.72, SE =0.67, P=0.01).

Island (Site
The MDS including all surveyed flowering stems showed a clear differentiation between islands with and without deer (stress =0.13), with an accompanying ANOSIM finding total floral communities varied more between islands grouped by deer presence than within island groups (R=0.627,P=0.01).The MDS including only 'usable' stems below 1.5 m demonstrated a greater difference in plant species community composition between islands with and without deer (stress=0.09,Figure 2 and Table S4) and again a significant ANOSIM result (R=0.702P=0.004).A complete list of all surveyed plants, their total stem counts and their native status in British Columbia can be found in the supplemental materials (Table S1).
Native 'usable' plant species richness from the entire field season was higher on unbrowsed islands (mean=10.5 ± 4.0) than browsed islands (mean=3.5±1.9,Table 1), and a Welch 2-sided t-test found them J o u r n a l P r e -p r o o f to significantly differ (t = 3.22, df = 3.97, p-value = 0.033).
Figure 2. MDS of 'usable' floral resources under 1.5 m placed on a 2-dimensional plane amongst islands with (yellow triangle) and without (blue square) deer present using the Bray-Curtis dissimilarity index.
Axes distances represent a visualization of differences in plant communities.Each different plant species is denoted with a black circle and the local status of the species as either Native (N) or Exotic (E).Drawings of plants highly associated with deer-absent (great camas, right) and deer-present islands (common dandelion, left) are featured on the figure.Stress = 0.09.

Bumblebee Abundances & Models
A total of 775 bees were surveyed (600 females and 175 males), with 88 kept as voucher specimens.A subset of 522 bumblebees observed actively foraging below 1.5 m were used for all analyses, including 6 different species (see Table 2 for female species breakdown).The excluded bees surveyed above 1.5 m were observed foraging mostly on tall ocean spray (Holodiscus discolor) and Himalayan blackberry (Rubus armeniacus) bushes and in the canopy of Arbutus (Arbutus menziesii) trees.Only queens were seen until May 5 ththe start of the second round of surveys and the day when the first worker was seen on Portland Islandafter which predominately workers were seen.The first males were seen on July 1 st on all islands surveyed that day.
The highest AIC ranked GLMM without uniformity or dispersion violations was selected to model total bumblebee abundance.For the female bumblebee analysis this was the second ranked model and included the log of mean stem count under 1.5 m, mean tree stand basal area, and mean temperature as fixed effect covariates and site ID as a random effect (Table S2).The GLMM found 'usable' stem count (P<0.0001) and temperature (P=0.07) to both have positive relationships with female bee abundances, with stem count significantly so.Mean tree stand basal area had a negative significant relationship with female abundances (P=0.03).The first ranked model for the male bumblebee GLMM was selected and included deer presence, mean temperature, and wind speed as fixed effects and site ID as a random effect.This model found deer presence and temperature both to have a significant negative relationship on male bee counts (P<0.000001, and P=0.0001 respectively) and wind speed to have a slight positive relationship with male bee counts (P=0.03)(Table S3).The PSEM confirmed that the presence of black-tailed deer has a direct negative relationship with abundance of flowering stems used by bumblebees, and an indirect negative relationship with female bumblebee abundances (independence claim P=0.54; Fisher's C=5.05).For this test the null hypothesis is that the model's data is different from the ideal fit, meaning a P>0.5 indicates a good representation of the data (Grace and Keeley, 2006).R 2 of values for both the female bee abundance GLMM and Linear Mixed Model were relatively high, indicating that no other important variables to be included were missed (R 2 Marginal =0.53, R 2 Conditional =0.55 and R 2 Marginal =0.26, R 2 Conditional =0.56 respectively, Figure 3).A PSEM was not conducted for male bumblebees as deer presence was a significant covariate and concurrent stem count was not found to be a significant covariate in the male GLMM.The MDS of female bumblebees shows that deer-absent islands have less variation in species composition amongst them, whereas deer-present islands have much greater variation (Figure 4, stress=0.08).The ANOSIM comparing bee assemblages between and within islands groups found there was no significant difference in community composition (R=0.19,P=0.1).Certain bee species had slight associations with deer absent or deer present islands (Figure 5), likely due to the plant species available on each island.Many native plant species were only found on deer absent islands, while other species (often exotic) were much more abundant on deer present islands (Table S1).For example, great camas (Camassia leichtlinii) was only found on deer absent islands, and the exotic common dandelion (Taraxacum officinale) was found at much greater abundances on islands with deer present.To understand plant-bee associations, the three most frequently foraged plants were identified for each bee species (Table 2).All bee species were most observed on native plants, except for B. vosnesenskii which was most often observed foraging from hairy cat's ear (Hypochaeris radicata) and ribwort plantain (Plantago lanceolata), exotic plant species closely associated with deer-present islands that were rarely visited by other bee species.Salal (Gaultheria shallon) and broadleaf stonecrop (Sedum spathulifolium) were the most utilized plants, with 4 of the 6 surveyed bee species foraging principally on either one (Table 2).B. sitkensis was the least observed species, with most sightings on Owl Island.The Western Bumblebee (B.occidentalis) was not observed on any sampled islands.
Table 2.All surveyed bumblebee species (foraging under 1.5 m) and the 3 plants they were found foraging on the most frequently.Percentage of foraging bees shows relative amount of total surveyed bees seen on the listed plant species.Only females were included in this count, as males were not identified to species.

Discussion
Our study contributes to the growing body of evidence on how hyperabundant herbivores can cause trophic cascade disruptions, altering native flora and fauna communities (Atkins et al., 2019;Louthan et al., 2019).We found that the presence of black-tailed deer on islands within the Salish Sea of British Columbia was associated with a decrease in flowering stems below 1.5 m (Figure 6), and that native plant species relied on by bumblebees were even more affected.We also found that bumblebee abundances were indirectly negatively impacted by hyperabundant deer populations via this reduction in floral resources.Floral resources that bees used were largely made up of native flowering herbaceous plants and shrubs, and were significantly reduced on islands with deer present.Previous research in the Salish Sea has documented negative impacts of deer hyperabundance on songbirds through the loss and modification of understory flora (Martin et al., 2011;Martin et al., 2014), and on herbaceous flora of cultural significance (Arcese et al., 2014;Gonzales & Arcese, 2008).Our results support the call for urgent deer management in order to protect and recover native plant and animal communities in the Salish Sea.J o u r n a l P r e -p r o o f We found deer density estimates to be much higher than previous studies in the Gulf Islands and believe this effect is most likely a result of preferred foraging habitat as opposed to growth in deer populations (Jung and Kukka, 2016).In Martin et al (2011) and Arcese et al (2014) deer estimates were calculated within forest habitats, whereas in our study we estimated deer density in open meadows which is a preferred habitat of deer (Collins and Urness, 1983), thus we would expect deer densities to be higher within these habitats.We acknowledge that these values do not represent actual deer densities for these entire islands, but instead site-specific densities which were adequate for comparison between sites for this study as all sites were located in similar habitat types.The high densities recorded in GOAE indicate the extreme browsing pressure these ecosystems are likely to be experiencing compared to other ecosystems on the same islands.Deer densities in predator-present locations across North America range from 0.03 to 0.3 deer/ha but averages around 0.12 (Ripple and Beschta, 2012), and align with recommended goals for density reduction in this area to reach 0.1 for ecological recovery (Martin et al., 2011).These results highlight the need for better methods for accurate density estimates to inform implementation of deer management and monitoring.

# of bees
Female bumblebee abundances were significantly and positively related to average stem count of 'usable' flowers in our GLMMs.'Usable' stem count is clearly an important covariate in this relationship as it was included in all compared female-bee GLMMs (which were all within 2 AIC of each other), with similar effect size throughout.The potential confounding variables island size and isolation were not included in most of the compared female-bee GLMMs, demonstrating they did not play a role in these results as might be expected by classic island biogeography theory.The PSEM showed that 'usable' stem counts had a direct relationship with observed female bumblebee abundances, and that deer presence had an indirect negative relationship with female bumblebees most likely due to their depletion of floral resources that bees rely on for survival and production.As bumblebees can travel distances over 1-2 km to forage from areas with available resources (Osborne et al., 2008), observations of foraging behaviour on each island does not necessarily indicate where the bumblebee's home nest site is located, but instead shows the relative importance of the resources in each study site to the general bumblebee population in the island archipelago (Hagen et al., 2011).
Male bumblebee counts were found to be significantly affected by deer presence in the GLMM, without concurrent stem count explaining this relationship.This may be due to male abundances being less impacted by the amount of current available flower stems, and more by the success of their colony that season, leading to greater male production (Pelletier and McNeil, 2003).Combining the results of both J o u r n a l P r e -p r o o f the male and female bumblebees, this study suggests that deer density may have short-term impacts on bumblebee communities through changes to active foraging and a longer-term impact through reduced reproductive success.However, for a better understanding of impacts to reproductive success, an analysis of newly produced queens should be included (Williams et al., 2012).The male GLMM results could also be in part due to the 'heat dome' during mid-June that occurred right before males began to emerge (Overland, 2021).The extreme heat dried up and killed many of the season's remaining flowering plants, leading to less of a difference of floral resources between islands with or without deer, which would explain why stem count was not significantly related to male abundances.Further study into the impacts of extreme weather on floral resources and pollinators are needed as these events are likely to continue with climate change (Overland, 2021).
Different bee species forage from slightly different plant assemblages based on their tongue length and resource requirements (Pellissier et al., 2013).The high use of native plants by the majority of bumblebee species demonstrates a crucial role for restoration and protection of remaining native flora throughout the area to support the bumblebee community.B. vosnesenskii was observed foraging in open meadows on exotic hairy cat's ear and ribwort plantain more than any other bumblebee species, potentially showing a greater resilience to floral community change caused by deer browsing than other species.B. vosnesenskii is a new arrival to the area, with its historic distribution further south in the United States (Fraser et al., 2012).Warming temperatures and the ability to forage on widespread exotic plant species may help explain its recent colonization of these Gulf islands.
It was hypothesised that B. occidentalis may be present at greater numbers on these small, remote islands compared to the mainland, as smaller, more isolated islands may present a refuge from invasive species and disease or parasite spread (Bennett and Arcese, 2013).Mite overloads have been identified as a major contributor to B. occidentalis declines, and are often spread from managed bees (Canadian Wildlife Service and Committee on the Status of Endangered Wildlife in Canada, 2014).Unfortunately, the redlisted B. occidentalis was not observed on any survey islands throughout the entire season.This demonstrates the decrease in B. occidentalis populations, which were historically one of the most common bumblebees in this area (Colla and Ratti, 2010).
Origins of plants were included in the floral MDS (Figure 2), showing that the majority of 'usable' plant species more closely associated with deer-present islands are exotic, while most plant species associated with deer-absent islands are native.While previous studies have shown that ungulates negatively impact bumblebees through browsing (Black et al., 2011;Filazzola et al., 2020;Louthan et al., 2019), few have investigated the overlap in preferred floral resources between deer and bees (DeBano et al., 2016).This study's results suggest that deer browsing is having a targeted impact on floral communities upon which bumblebees rely.
While attempts were made to select a range of islands with and without deer present, some differences were unavoidable.D'Arcy Island and Portland Island arguably have the most differing histories in regards to deer presence, with D'Arcy Island having a more recent known deer colonization in the 1960s rather than a long history of deer always being present.Portland, on the other hand, represents the opposite situation, where deer and sheep inhabited the island until the 1990s when they were removed (personal correspondence, P. Linton).While passive restoration of the plant communities on Portland are underway, much of the Garry oak ecosystem there resembles browsed communities, with a lower than average floral stem count than other islands without deer present and greater presence of invasive species.
Our results demonstrate that deer have had a clear impact on floral communities used by bumblebees.Field observations also show an unappreciated usefulness of certain invasive species such as Scotch broom (Cytisus scoparius) for pollinators.For example, because Scotch broom is resistant to deer browsing it provides floral resources for bumblebees in otherwise resource-poor areas, while also J o u r n a l P r e -p r o o f providing refuge for more fragile plant species in its undergrowth that have been completely browsed by deer elsewhere (Simberloff et al., 2003).Many of these islands are under the protection of Parks Canada, where there is currently minimal active management, with broom removal being one of the only programs run on these smaller islands (Parks Canada, 2006).While there are negative impacts from the spread of this invasive species including reduced soil phosphorous, increased fire risks, and suppression of native flora (Carter et al., 2019;Shaben and Myers, 2010), our results show that this management program may be placing bumblebees under even greater stress for finding floral resources on islands where deer remain.
Restoring these ecosystems may only be effective if broom removal is done in accompaniment with deer reductions and active restoration of native species, otherwise very few native plant species will have a chance of re-establishing to support pollinators.
Restoration of these areas is slow and difficult (Shackelford et al., 2019).Without long-term, sustained efforts, areas absent of deer will not easily return to previous meadows filled with all the native flora that were once so prevalent (Clements, 2013).Although restoration of the ecological integrity of the area would require a substantive investment of resources, bumblebee populations would likely respond quickly and positively.This study's results shows deer-absent islands are currently acting as areas of refuge for bumblebee populations; due to bumblebee's short colony life-span of one year, queens are establishing new nests every year allowing continual colonization to deer populated islands (Koch et al., 2012).However, with continued un-checked deer populations, Garry oak ecosystems and their support of bumblebee populations are at risk.Current methods of placing GOAEs in protected areas is not enough these areas must also be actively managed to reduce browsing pressure (Martin et al., 2011).With over-browsing from hyperabundant deer many plants that bees depend on have no opportunity to flower, set seed, or reproduce (Beschta and Ripple, 2009).Endangered plant communities and wild pollinators often rely on each other, where the health and persistence of one depends on the other (Koh et al., 2004;Vieira et al., 2013).With continued bumblebee declines along with widespread habitat reduction and degradation, it is imperative to properly protect their remaining habitat for their persistence while populations are still able to recover.

Conclusion
Hyperabundant black-tailed deer are causing significant negative impacts on the Garry Oak and Associated Ecosystems of British Columbia.Through over-browsing, black-tailed deer are causing a complex cascading reaction through the ecological community and an overall shift in ecosystem composition, with specific ramifications for native plant species and bumblebee populations.These results show that if restoration of diverse native plant communities and their associated pollinators in GOAE is a conservation goal, then determining the most cost-effective deer management options that meet both ecological and societal concerns is a research priority.

Figure 1 .
Figure 1.Map of the 10 islands with study sites.Surveys were conducted on islands outlined in black.Survey sites with deer present are indicated by filled dots (browsed), and sites without deer present are indicated by empty dots (unbrowsed).
r n a l P r e -p r o o f

Figure 3 .
Figure 3. Piece-wise structural equation modelling results.Letters (a) through (d) show model estimates and standard errors, with the accompanying table showing relative effects of each.Red arrows represent negative relationships and black arrows represent positive relationships.R 2 M stands for Marginal R 2 and R 2 C stands for Conditional R 2 .P = 0.537, which in this case indicates the model represents the data well.
18% increase in mean bumblebee count per 1 degree increase in temperature (°C) (b) 245% decrease in the log of mean 'usable' stem count (per m 2 ) on deer-present islands (c) 6% decrease in mean bumblebee count per 1 additional tree in basal area count (d) 68% increase in mean bumblebee count per 1 additional unit of log(mean 'usable' stems) Relative Effect J o u r n a l P r e -p r o o f

Figure 6 .
Figure 6.Surveying floral resources on islands with (left) and without (right) black-tailed deer in the spring of 2021.Photos taken on the same day (April 19) within the same ecosystem type.