Final report on bee biodiversity across Denver Parks

Author

John M. Mola, Nicki Bailey, Nathan Comai, Nancy Bartholomew

Purpose

The purpose of this report is to provide a “heat map” of different dimensions of bee diversity across Denver Parks. Additionally, we provide tables of interest related to bee richness, abundance, and diversity across parks as well as plants that bees were commonly captured from.

Future analyses are expected but are beyond the scope of the initial proposal. However, these are likely to be of great interest to Denver Parks and Recreation (DPR) and will be reported to DPR when ready. These include:

  • Analysis of biodiversity with respect to local (e.g. park typology, plant diversity, honeybee abundance, etc) and landscape scale factors (e.g. impervious surface, residential density, income, etc)
  • Detailed plant analysis
  • Guild-specific pollinator studies (e.g. proportion of native/non-native, parasitic, social, etc)
  • Butterfly abundance, richness, and diversity data

Please let us know if there are any other visualizations or analyses of immediate interest.

Show the code
# WRANGLING TO FIND RICHNESS ----------------------------------------------

# total count of each species across parks
df_bee_count <- df_full_bees %>% 
  group_by(park_name, id_code, genus, species) %>% 
  summarise(n_individuals = n())

# summarized count of bee species richness across parks
df_sum_bee_count <- df_bee_count %>% 
  group_by(park_name, id_code) %>% 
  summarise(richness = n())

# calculating basic diversity indices on bee data for each site; this does not contain the park name, but does include the location code
df_diversity <- df_bee_count %>% 
  # these steps needed for vegan's format
  mutate(species = paste(genus, species)) %>% 
  select(-genus) %>% 
  ungroup() %>% 
  group_by(id_code) %>% 
  summarise(total_count = sum(n_individuals), 
            shannon.di = diversity(n_individuals, index = "shannon"), 
            simpson.di = diversity(n_individuals, index = "simpson")) %>% 
  arrange(desc(shannon.di))


# JOIN BEE RICHNESS TO PARK SPATIAL DATA ----------------------------------

# joins the summarized counts to the shapefile

df_shp_join <- df_shp_parks %>% 
  inner_join(., df_sum_bee_count, by = c("LOC_CODE" = "id_code")) %>% 
  inner_join(., df_diversity, by = c("LOC_CODE" = "id_code"))

Map of Species Richness (i.e. number of bee species observed) Across Parks

Show the code
p_map_richness<- ggplot() +
  geom_sf(data = shp_nbr, fill = "grey90", color = "grey70") +
  # geom_sf(data = shp_hud, aes(fill = LOWMODPCT)) +
  #geom_sf(data = shp_tree, aes(fill = PCT_CANOPY)) +
  
  #new_scale_fill()+
  # geom_sf(data = shp_flt_parks, aes(color = PARK_CLASS)) +
  # geom_sf(data = df_shp_join, aes(fill = richness, text = park_name)) +
  geom_point(data = df_shp_join, aes(color = richness, text = park_name, x = LONGITUDE, y = LATITUDE), size = 3, alpha = 0.7) +
  scale_color_viridis_c() +
  #labs(fill = "Overall Neighborhood Equity Score", x = "Longitude", y = "Latitude", color = "") +
  theme_bw(base_family = "nunito") +
  theme()

ggplotly(p_map_richness)

Map of Bee Abundance Across Parks

Show the code
p_map_abundance<- ggplot() +
  geom_sf(data = shp_nbr, fill = "grey90", color = "grey70") +
  # geom_sf(data = shp_hud, aes(fill = LOWMODPCT)) +
  #geom_sf(data = shp_tree, aes(fill = PCT_CANOPY)) +
  
  #new_scale_fill()+
  # geom_sf(data = shp_flt_parks, aes(color = PARK_CLASS)) +
  # geom_sf(data = df_shp_join, aes(fill = richness, text = park_name)) +
  geom_point(data = df_shp_join, aes(color = total_count, text = park_name, x = LONGITUDE, y = LATITUDE), size = 3, alpha = 0.7) +
  scale_color_viridis_c() +
  #labs(fill = "Overall Neighborhood Equity Score", x = "Longitude", y = "Latitude", color = "") +
  theme_bw(base_family = "nunito") +
  theme()

ggplotly(p_map_abundance)

Bee richness vs. Bee abundance

City Park and Platte Farm Open space have perhaps more species than we might expect given the number of specimens collected. Washington Park and Argo Park may have fewer species than expected given the number of individuals observed.

Show the code
p_abunrich <- ggplot(data = df_shp_join, aes(x = total_count, y = richness, text = park_name)) +
  geom_point(size = 3) +
  # geom_abline(intercept = 0, slope = 1) +
  # geom_smooth(method = "lm") +
  theme_classic() +
  labs(x = "Number of Specimens", y = "Species Richness")

ggplotly(p_abunrich)

Map of Species Diversity Across Parks

Using Shannon Diversity Index; helps us understand richness with some respect to individual abundance. Higher numbers = more diverse. (This is strongly correlated with other diversity indices in our data)

Show the code
p_map_diversity <- ggplot() +
  geom_sf(data = shp_nbr, fill = "grey90", color = "grey70") +
  # geom_sf(data = shp_hud, aes(fill = LOWMODPCT)) +
  #geom_sf(data = shp_tree, aes(fill = PCT_CANOPY)) +
  
  #new_scale_fill()+
  # geom_sf(data = shp_flt_parks, aes(color = PARK_CLASS)) +
  # geom_sf(data = df_shp_join, aes(fill = richness, text = park_name)) +
  geom_point(data = df_shp_join, aes(color = shannon.di, text = park_name, x = LONGITUDE, y = LATITUDE), size = 3, alpha = 0.7) +
  scale_color_viridis_c() +
  #labs(fill = "Overall Neighborhood Equity Score", x = "Longitude", y = "Latitude", color = "") +
  theme_bw(base_family = "nunito") +
  theme()

ggplotly(p_map_diversity)

Summary Table of Richness, Abundance, and Diversity Across Parks

Show the code
df_plant_use <- df_full_bees %>% 
  filter(!is.na(floral_association_genspe)) %>%   
  mutate(floral_species = paste(floral_association_genus, floral_association_species)) %>%  
  group_by(id_code, floral_species) %>% 
  count() %>% 
  group_by(id_code) %>% 
  count(name = "plant_use")

df_shp_join %>% 
  inner_join(., df_plant_use, by = c("LOC_CODE" = "id_code")) %>% 
  st_drop_geometry() %>% 
  dplyr::select(park_name, richness, total_count, shannon.di, simpson.di, plant_use) %>% 
  mutate_if(is.numeric, ~round(., 2)) %>% 
  arrange(park_name) %>% 
  DT::datatable(colnames = c("Park Name", "Richness", "Number of Individuals", "Shannon Diversity", "Simpson Diversity", "Plants used by bees"), 
                extensions = "Buttons", 
                options = list(iDisplayLength = 25, 
                               dom = 'Bfrtip',
                               buttons = c('copy', 'csv', 'excel', 'pdf')))

Top Plants Used by Bees (Lumped across all parks)

Some names may not be merged or taxonomy may be uncorrected if species names have changed. Data is preliminary.

Show the code
df_full_bees %>% 
  filter(!is.na(floral_association_genspe)) %>% 
  mutate(floral_species = paste(floral_association_genus, floral_association_species)) %>% 
  count(floral_species) %>% 
  arrange(desc(n)) %>% 
  DT::datatable(colnames = c("Floral Species", "Number of bees collected from that plant"), 
                extensions = "Buttons", 
                options = list(iDisplayLength = 20, 
                               dom = 'Bfrtip',
                               buttons = c('copy', 'csv', 'excel', 'pdf')))

Plants used by bees, all species, all parks

Some names may not be merged or taxonomy may be uncorrected if species names have changed. Data is preliminary.

Show the code
df_full_bees %>% 
  filter(!is.na(floral_association_genspe)) %>%   
  mutate(floral_species = paste(floral_association_genus, floral_association_species)) %>%  
  group_by(park_name, floral_species) %>% 
  count() %>% 
  DT::datatable(colnames = c("Park", "Floral Species", "Number of individual bees captured on plant"), 
                extensions = "Buttons", 
                options = list(iDisplayLength = 20, 
                               dom = 'Bfrtip',
                               buttons = c('copy', 'csv', 'excel', 'pdf')))