# R script for "Nutritional composition of pollen stores in managed bees across 
# European agro-ecosystems reveals species-specific differences but 
# limited pesticide effects"

# Gekière et al. (2026)

# Published in Ecological Entomology

setwd("XXX")

library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggeffects)
library(lme4)
library(car)
library(glmmTMB)
library(DHARMa)
library(emmeans)
library(multcomp)
library(ggpubr)

######### LOAD DATASET #########

data <- read_excel("XXX")
data$CountryID <- as.factor(data$CountryID)
data$SiteID <- as.factor(data$SiteID)
data$Crop <- as.factor(data$Crop)
data$Type <- as.factor(data$Type)
data$fl_prop_mfl <- 1 - data$fl_prop_nmfl ## Proportion of pollen from crops
data$logTU <- log(data$MixtureTU_avg_pollen.stores + 0.1) ## Take log of pesticide risk

# Subset dataset per bee species
data_apis <- subset(data, Type == "apis")
data_bombus <- subset(data, Type == "bombus")
data_osmia <- subset(data, Type == "osmia")

######### INTERSP SUMMARY #########

data %>% group_by(Type) %>%
  summarise(n = n(),
            mean_prot = mean(protein, na.rm = T),
            mean_lipid = mean(lipid, na.rm = T),
            mean_PL = mean(p2l, na.rm = T))

######### PROTEIN ANALYSES #########

# Model
model_prot <- lmer(data = data, protein ~ Type*logTU + Crop + (1|CountryID/SiteID))
Anova(model_prot, type = "II")
simulationOutput <- simulateResiduals(fittedModel = model_prot, plot = T)

# Emmeans to test differences between species
(emm <- emmeans(model_prot, pairwise ~ Type, adjust = "fdr"))
cld(emm, alpha = 0.05, Letters = LETTERS, adjust = "fdr")

# Emtrends to test effect of logTU per species
test(emtrends(model_prot, pairwise ~ Type, var = "logTU", adjust = "fdr"))

# Extract model prediction for comparison among species
pred_inter_prot <- ggemmeans(model_prot, terms = "Type", type = "fixed")

# Plot for comparison among species
(inter_graph_prot <- ggplot(data, aes(x = Type, y = protein)) +
    geom_jitter(aes(fill = CountryID, shape = Crop),
                width = 0.2, alpha = 0.5, size = 5, colour = "black") +
    geom_point(data = pred_inter_prot, aes(x = x, y = predicted),
               size = 7, fill = "black", inherit.aes = FALSE) +
    geom_errorbar(data = pred_inter_prot,
                  aes(x = x, ymin = conf.low, ymax = conf.high),
                  width = 0, size = 3, inherit.aes = FALSE) +
    labs(x = "Bee species", y = "Protein content (µg / mg)") +
    ylim(0, 570) +
    scale_x_discrete(labels = c("Apis", "Bombus", "Osmia")) +
    scale_shape_manual(values = c(21,24)) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size = 20),
          axis.text.x = element_text(size = 16, face = "italic"),
          axis.text.y = element_text(size = 16),
          plot.tag = element_text(size = 25, face = "bold"),
          legend.position = "none",
          plot.margin = margin(5, 5, 5, 5)))

# Extract model prediction across logTU for Apis
pred_apis_prot <- ggemmeans(model_prot, terms = c("logTU [-2.13:8 by=0.1]", "Type [apis]"), 
                            type = "fixed")

# Plot across logTU for Apis
(apis_graph_prot <- ggplot(data = data_apis, aes(x = logTU, y = protein)) +
    geom_point(aes(fill = CountryID, shape = Crop), alpha = 0.5, size = 5, colour = "black") +
    geom_line(data = pred_apis_prot, aes(x = x, y = predicted), 
              colour = "black", size = 4) +
    geom_ribbon(data = pred_apis_prot, aes(x = x, ymin = conf.low, ymax = conf.high), 
                inherit.aes = FALSE, alpha = 0.2, fill = "grey20") +
    labs(x = "Pesticide risk", y = "Protein content (µg / mg)") +
    scale_shape_manual(values = c(21,24)) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size = 20),
          axis.text.x = element_text(size = 16),
          axis.text.y = element_text(size = 16),
          plot.tag = element_text(size = 25, face = "bold"),
          legend.position = "none"))

# Extract model prediction across logTU for Bombus
pred_bombus_prot <- ggemmeans(model_prot, terms = c("logTU [-2.31:11.38 by=0.1]", "Type [bombus]"), 
                            type = "fixed")

# Plot across logTU for Bombus
(bombus_graph_prot <- ggplot(data = data_bombus, aes(x = logTU, y = protein)) +
    geom_point(aes(fill = CountryID, shape = Crop), alpha = 0.5, size = 5, colour = "black") +
    geom_line(data = pred_bombus_prot, aes(x = x, y = predicted), 
              colour = "black", size = 4) +
    geom_ribbon(data = pred_bombus_prot, aes(x = x, ymin = conf.low, ymax = conf.high), 
                inherit.aes = FALSE, alpha = 0.2, fill = "grey20") +
    labs(x = "Pesticide risk", y = "Protein content (µg / mg)") +
    scale_shape_manual(values = c(21,24)) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size = 20),
          axis.text.x = element_text(size = 16),
          axis.text.y = element_text(size = 16),
          plot.tag = element_text(size = 25, face = "bold"),
          legend.position = "none"))

# Extract model prediction across logTU for Osmia
pred_osmia_prot <- ggemmeans(model_prot, terms = c("logTU [-2.31:7.6 by=0.1]", "Type [osmia]"), 
                              type = "fixed")

# Plot across logTU for Osmia
(osmia_graph_prot <- ggplot(data = data_osmia, aes(x = logTU, y = protein)) +
    geom_point(aes(fill = CountryID, shape = Crop), alpha = 0.5, size = 5, colour = "black") +
    geom_line(data = pred_osmia_prot, aes(x = x, y = predicted), 
              colour = "black", size = 4) +
    geom_ribbon(data = pred_osmia_prot, aes(x = x, ymin = conf.low, ymax = conf.high), 
                inherit.aes = FALSE, alpha = 0.2, fill = "grey20") +
    labs(x = "Pesticide risk", y = "Protein content (µg / mg)") +
    scale_shape_manual(values = c(21,24)) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size = 20),
          axis.text.x = element_text(size = 16),
          axis.text.y = element_text(size = 16),
          plot.tag = element_text(size = 25, face = "bold"),
          legend.position = "none"))

######### LIPID ANALYSES #########

# Model
model_lip <- lmer(data = data, lipid ~ Type*logTU + Crop + (1|CountryID/SiteID))
Anova(model_lip, type = "II")
simulationOutput <- simulateResiduals(fittedModel = model_lip, plot = T)

# Emmeans to test differences between species
(emm <- emmeans(model_lip, pairwise ~ Type, adjust = "fdr"))
cld(emm, alpha = 0.05, Letters = LETTERS, adjust = "fdr")

# Emtrends to test effect of logTU per species
test(emtrends(model_lip, pairwise ~ Type, var = "logTU", adjust = "fdr"))

# Extract model prediction for comparison among species
pred_inter_lip <- ggemmeans(model_lip, terms = "Type", type = "fixed")

# Plot for comparison among species
(inter_graph_lip <- ggplot(data, aes(x = Type, y = lipid)) +
    geom_jitter(aes(fill = CountryID, shape = Crop),
                width = 0.2, alpha = 0.5, size = 5, colour = "black") +
    geom_point(data = pred_inter_lip, aes(x = x, y = predicted),
               size = 7, fill = "black", inherit.aes = FALSE) +
    geom_errorbar(data = pred_inter_lip,
                  aes(x = x, ymin = conf.low, ymax = conf.high),
                  width = 0, size = 3, inherit.aes = FALSE) +
    labs(x = "Bee species", y = "Lipid content (µg / mg)") +
    ylim(0, 150) +
    scale_x_discrete(labels = c("Apis", "Bombus", "Osmia")) +
    scale_shape_manual(values = c(21,24)) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size = 20),
          axis.text.x = element_text(size = 16, face = "italic"),
          axis.text.y = element_text(size = 16),
          plot.tag = element_text(size = 25, face = "bold"),
          legend.position = "none",
          plot.margin = margin(5, 5, 5, 5)))

# Extract model prediction across logTU for Apis
pred_apis_lip <- ggemmeans(model_lip, terms = c("logTU [-2.13:8 by=0.1]", "Type [apis]"), 
                            type = "fixed")

# Plot across logTU for Apis
(apis_graph_lip <- ggplot(data = data_apis, aes(x = logTU, y = lipid)) +
    geom_point(aes(fill = CountryID, shape = Crop), alpha = 0.5, size = 5, colour = "black") +
    geom_line(data = pred_apis_lip, aes(x = x, y = predicted), 
              colour = "black", size = 4) +
    geom_ribbon(data = pred_apis_lip, aes(x = x, ymin = conf.low, ymax = conf.high), 
                inherit.aes = FALSE, alpha = 0.2, fill = "grey20") +
    labs(x = "Pesticide risk", y = "Lipid content (µg / mg)") +
    scale_shape_manual(values = c(21,24)) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size = 20),
          axis.text.x = element_text(size = 16),
          axis.text.y = element_text(size = 16),
          plot.tag = element_text(size = 25, face = "bold"),
          legend.position = "none"))

# Extract model prediction across logTU for Bombus
pred_bombus_lip <- ggemmeans(model_lip, terms = c("logTU [-2.31:11.38 by=0.1]", "Type [bombus]"), 
                              type = "fixed")

# Plot across logTU for Bombus
(bombus_graph_lip <- ggplot(data = data_bombus, aes(x = logTU, y = lipid)) +
    geom_point(aes(fill = CountryID, shape = Crop), alpha = 0.5, size = 5, colour = "black") +
    geom_line(data = pred_bombus_lip, aes(x = x, y = predicted), 
              colour = "black", size = 4) +
    geom_ribbon(data = pred_bombus_lip, aes(x = x, ymin = conf.low, ymax = conf.high), 
                inherit.aes = FALSE, alpha = 0.2, fill = "grey20") +
    labs(x = "Pesticide risk", y = "Lipid content (µg / mg)") +
    scale_shape_manual(values = c(21,24)) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size = 20),
          axis.text.x = element_text(size = 16),
          axis.text.y = element_text(size = 16),
          plot.tag = element_text(size = 25, face = "bold"),
          legend.position = "none"))

# Extract model prediction across logTU for Osmia
pred_osmia_lip <- ggemmeans(model_lip, terms = c("logTU [-2.31:7.6 by=0.1]", "Type [osmia]"), 
                             type = "fixed")

# Plot across logTU for Osmia
(osmia_graph_lip <- ggplot(data = data_osmia, aes(x = logTU, y = lipid)) +
    geom_point(aes(fill = CountryID, shape = Crop), alpha = 0.5, size = 5, colour = "black") +
    geom_line(data = pred_osmia_lip, aes(x = x, y = predicted), 
              colour = "black", size = 4) +
    geom_ribbon(data = pred_osmia_lip, aes(x = x, ymin = conf.low, ymax = conf.high), 
                inherit.aes = FALSE, alpha = 0.2, fill = "grey20") +
    labs(x = "Pesticide risk", y = "Lipid content (µg / mg)") +
    scale_shape_manual(values = c(21,24)) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size = 20),
          axis.text.x = element_text(size = 16),
          axis.text.y = element_text(size = 16),
          plot.tag = element_text(size = 25, face = "bold"),
          legend.position = "none"))

######### PROTEIN-TO-LIPID ANALYSES #########

# Model
model_PL <- lmer(data = data, p2l ~ Type*logTU + Crop + (1|CountryID/SiteID))
Anova(model_PL, type = "II")
simulationOutput <- simulateResiduals(fittedModel = model_PL, plot = T)

# Emmeans to test differences between species
(emm <- emmeans(model_PL, pairwise ~ Type, adjust = "fdr"))
cld(emm, alpha = 0.05, Letters = LETTERS, adjust = "fdr")

# Emtrends to test effect of logTU per species
test(emtrends(model_PL, pairwise ~ Type, var = "logTU", adjust = "fdr"))

# Extract model prediction for comparison among species
pred_inter_PL <- ggemmeans(model_PL, terms = "Type", type = "fixed")

# Plot for comparison among species
(inter_graph_PL <- ggplot(data, aes(x = Type, y = p2l)) +
    geom_jitter(aes(fill = CountryID, shape = Crop),
                width = 0.2, alpha = 0.5, size = 5, colour = "black") +
    geom_point(data = pred_inter_PL, aes(x = x, y = predicted),
               size = 7, fill = "black", inherit.aes = FALSE) +
    geom_errorbar(data = pred_inter_PL,
                  aes(x = x, ymin = conf.low, ymax = conf.high),
                  width = 0, size = 3, inherit.aes = FALSE) +
    labs(x = "Bee species", y = "P:L ratio") +
    ylim(0, 25) +
    scale_x_discrete(labels = c("Apis", "Bombus", "Osmia")) +
    scale_shape_manual(values = c(21,24)) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size = 20),
          axis.text.x = element_text(size = 16, face = "italic"),
          axis.text.y = element_text(size = 16),
          plot.tag = element_text(size = 25, face = "bold"),
          legend.position = "none",
          plot.margin = margin(5, 5, 5, 5)))

# Extract model prediction across logTU for Apis
pred_apis_PL <- ggemmeans(model_PL, terms = c("logTU [-2.13:8 by=0.1]", "Type [apis]"), 
                           type = "fixed")

# Plot across logTU for Apis
(apis_graph_PL <- ggplot(data = data_apis, aes(x = logTU, y = p2l)) +
    geom_point(aes(fill = CountryID, shape = Crop), alpha = 0.5, size = 5, colour = "black") +
    geom_line(data = pred_apis_PL, aes(x = x, y = predicted), 
              colour = "black", size = 4) +
    geom_ribbon(data = pred_apis_PL, aes(x = x, ymin = conf.low, ymax = conf.high), 
                inherit.aes = FALSE, alpha = 0.2, fill = "grey20") +
    labs(x = "Pesticide risk", y = "P:L ratio") +
    scale_shape_manual(values = c(21,24)) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size = 20),
          axis.text.x = element_text(size = 16),
          axis.text.y = element_text(size = 16),
          plot.tag = element_text(size = 25, face = "bold"),
          legend.position = "none"))

# Extract model prediction across logTU for Bombus
pred_bombus_PL <- ggemmeans(model_PL, terms = c("logTU [-2.31:11.38 by=0.1]", "Type [bombus]"), 
                             type = "fixed")

# Plot across logTU for Bombus
(bombus_graph_PL <- ggplot(data = data_bombus, aes(x = logTU, y = p2l)) +
    geom_point(aes(fill = CountryID, shape = Crop), alpha = 0.5, size = 5, colour = "black") +
    geom_line(data = pred_bombus_PL, aes(x = x, y = predicted), 
              colour = "black", size = 4) +
    geom_ribbon(data = pred_bombus_PL, aes(x = x, ymin = conf.low, ymax = conf.high), 
                inherit.aes = FALSE, alpha = 0.2, fill = "grey20") +
    labs(x = "Pesticide risk", y = "P:L ratio") +
    scale_shape_manual(values = c(21,24)) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size = 20),
          axis.text.x = element_text(size = 16),
          axis.text.y = element_text(size = 16),
          plot.tag = element_text(size = 25, face = "bold"),
          legend.position = "none"))

# Extract model prediction across logTU for Osmia
pred_osmia_PL <- ggemmeans(model_PL, terms = c("logTU [-2.31:7.6 by=0.1]", "Type [osmia]"), 
                            type = "fixed")

# Plot across logTU for Osmia
(osmia_graph_PL <- ggplot(data = data_osmia, aes(x = logTU, y = p2l)) +
    geom_point(aes(fill = CountryID, shape = Crop), alpha = 0.5, size = 5, colour = "black") +
    geom_line(data = pred_osmia_PL, aes(x = x, y = predicted), 
              colour = "black", size = 4) +
    geom_ribbon(data = pred_osmia_PL, aes(x = x, ymin = conf.low, ymax = conf.high), 
                inherit.aes = FALSE, alpha = 0.2, fill = "grey20") +
    labs(x = "Pesticide risk", y = "P:L ratio") +
    scale_shape_manual(values = c(21,24)) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size = 20),
          axis.text.x = element_text(size = 16),
          axis.text.y = element_text(size = 16),
          plot.tag = element_text(size = 25, face = "bold"),
          legend.position = "none"))

######### CHECK PESTICIDE RISK VS PROPORTION CROP IN POLLEN #########

ggplot(data = data, aes(x = fl_prop_mfl, y = logTU)) +
  geom_point() + geom_smooth(method = "lm")
cor.test(data$fl_prop_mfl, data$logTU, method = "pearson", use = "complete.obs")

######### COMBINE PLOTS FOR COMPARISON AMONG SPECIES (FIG 2) #########

ggarrange(inter_graph_prot, inter_graph_lip, inter_graph_PL,
          nrow = 3, align = "hv")

######### COMBINE PLOTS FOR PROT AND LIP WITHIN EACH SPECIES (FIG 3) #########

ggarrange(ggarrange(apis_graph_prot, apis_graph_lip, nrow = 1),
          ggarrange(bombus_graph_prot, bombus_graph_lip, nrow = 1),
          ggarrange(osmia_graph_prot, osmia_graph_lip, nrow = 1),
          nrow = 3, align = "hv")

######### COMBINE PLOTS FOR PL RATIO WITHIN EACH SPECIES (FIG 4) #########

ggarrange(apis_graph_PL, bombus_graph_PL, osmia_graph_PL,
          nrow = 3, align = "hv")