#Code to run analyses used in Climate-driven phenological shifts in emergence dates of British bees.
#Author: Chris Wyver
#Last Updated: 21/6/2023

####Load Libraries####
library(readr)
library(nlme)
library(performance)
library(dplyr)
library(tidyverse)
library(purrr)
library(tidyr)
library(broom)
library(FSA)
library(ggpubr)
library(cowplot)

####Load Data####
Emergence <- read_csv("Wyver_et_al_2023_Dataset/EmergenceDates.csv")
Traits <- subset(Emergence, select = c(1, 6:10))
Traits <- unique(Traits)

####Year Models####
yearfirst <- lme(Emergence ~ year + n + Northing, random= ~ 1|Species, data=Emergence)
summary(yearfirst)
r2_nakagawa(yearfirst, tolerance = 1e-330)

models <- Emergence %>% 
  nest(data = -Species) %>% 
  mutate(
    model = purrr::map(data, ~ lm(Emergence ~ year + n + Northing, data = .)),
    tidied = purrr::map(model, tidy),
    r_squared = map_dbl(model, ~ summary(.)$r.squared), 
    p_value_Year = map_dbl(model, ~ tidy(summary(.)) %>% filter(term == "year") %>% pull(p.value)),
    p_value_n = map_dbl(model, ~ tidy(summary(.)) %>% filter(term == "n") %>% pull(p.value)),
    p_value_Northing = map_dbl(model, ~ tidy(summary(.)) %>% filter(term == "Northing") %>% pull(p.value)),
    corrected_p_value_Year = p.adjust(p_value_Year, method = "BH"),
    corrected_p_value_n = p.adjust(p_value_n, method = "BH"),
    corrected_p_value_Northing = p.adjust(p_value_Northing, method = "BH")
  )

# View the tidied model output
tidy_models <- models %>% 
  select(Species, tidied, r_squared, corrected_p_value_Year, corrected_p_value_n, corrected_p_value_Northing) %>% 
  unnest(tidied)

tidy_models


#Get year stats
tidy_models_year <- subset(tidy_models, term == "year")
tidy_models_year_sig <- subset(tidy_models_year, corrected_p_value_Year < 0.05)

#Get n stats
tidy_models_n <- subset(tidy_models, term == "n")
tidy_models_n_sig <- subset(tidy_models_n, corrected_p_value_n < 0.05)

#Get latitude stats
tidy_models_Northing <- subset(tidy_models, term == "Northing")
tidy_models_Northing_sig <- subset(tidy_models_Northing, corrected_p_value_Northing < 0.05)

#Plot Figure 1
tidy_models_year<-merge(tidy_models_year, Traits, by = c("Species"))
first_models_plot<-tidy_models_year %>% mutate(Species = fct_reorder(Species, estimate))
first_models_plot$Genus <- word(first_models_plot$Species,1)
ggplot(first_models_plot)+
  geom_point(aes(x = estimate, y = Species, color = Genus))+
  geom_segment(aes(x = estimate-std.error, y = Species, xend = estimate+std.error, yend = Species, color = Genus))+
  theme_bw()+
  theme(axis.text=element_text(size=7),axis.title=element_text(size=15), axis.text.y = element_text(face = 'italic'), legend.text = element_text(face = 'italic'))+
  labs(x = "Change in emergence date (Days/Year)")+
  geom_vline(xintercept = 0)

#Plot Figure S1
ggplot(Emergence,aes(x = year, y = Emergence))+
  geom_point()+
  geom_smooth(method = "lm")+
  facet_wrap(~Species)+
  theme_bw()+
  theme(axis.text = element_text(size = 5))+
  theme(strip.text = element_text(face = "italic", size = 6.5))+
  labs(x = "Year", y = "Emergence Date (Day of Year)")



####Climate Models####
TGfirst <- lme(Emergence ~ TG + n + Northing, random=~1|Species, data=Emergence, na.action = na.omit)
summary(TGfirst)
r2_nakagawa(TGfirst, tolerance = 1e-330)

TGmodels <- Emergence %>% 
  nest(data = -Species) %>% 
  mutate(model = purrr::map(data, ~ lm(Emergence ~ TG + n + Northing, data = .)),
         tidied = purrr::map(model, tidy),
         r_squared = map_dbl(model, ~ summary(.)$r.squared), 
         p_value_TG = map_dbl(model, ~ tidy(summary(.)) %>% filter(term == "TG") %>% pull(p.value)),
         p_value_n = map_dbl(model, ~ tidy(summary(.)) %>% filter(term == "n") %>% pull(p.value)),
         p_value_Northing = map_dbl(model, ~ tidy(summary(.)) %>% filter(term == "Northing") %>% pull(p.value)),
         corrected_p_value_TG = p.adjust(p_value_TG, method = "BH"),
         corrected_p_value_n = p.adjust(p_value_n, method = "BH"),
         corrected_p_value_Northing = p.adjust(p_value_Northing, method = "BH"))

# View the tidied model output
TG_tidy_models <- TGmodels %>% 
  select(Species, tidied, r_squared, corrected_p_value_TG, corrected_p_value_n, corrected_p_value_Northing) %>% 
  unnest(tidied)

#Get year stats
tidy_models_TG <- subset(TG_tidy_models, term == "TG")
tidy_models_TG_sig <- subset(tidy_models_TG, corrected_p_value_TG < 0.05)

#Get n stats
tidy_models_n_TG <- subset(TG_tidy_models, term == "n")
tidy_models_n_TG_sig <- subset(tidy_models_n_TG, corrected_p_value_n < 0.05)

#Get latitude stats
tidy_models_Northing_TG <- subset(TG_tidy_models, term == "Northing")
tidy_models_Northing_TG_sig <- subset(tidy_models_Northing_TG, corrected_p_value_Northing < 0.05)

#Plot Figure 2
TG_models_plot<-tidy_models_TG %>% mutate(Species = fct_reorder(Species, estimate))
TG_models_plot$Genus <- word(TG_models_plot$Species,1)
ggplot(TG_models_plot)+
  geom_point(aes(x = estimate, y = Species, color = Genus))+
  geom_segment(aes(x = estimate-std.error, y = Species, xend = estimate+std.error, yend = Species, color = Genus))+
  theme_bw()+
  theme(axis.text=element_text(size=7),axis.title=element_text(size=15), axis.text.y = element_text(face = 'italic'), legend.text = element_text(face = 'italic'))+
  labs(x = "Change in emergence date (Days/°C)")+
  geom_vline(xintercept = 0)

#Plot Figure S2
ggplot(Emergence,aes(x = TG, y = Emergence))+
  geom_point()+
  geom_smooth(method = "lm")+
  facet_wrap(~Species)+
  theme_bw()+
  theme(axis.text = element_text(size = 5))+
  theme(strip.text = element_text(face = "italic", size = 6.5))+
  labs(x = "Temperature during 90-day climate window (°C)", y = "Emergence Date (Day of Year)")

####Traits Comparison####
Traits %>%
  group_by(Trait.Group) %>%
  summarise(count=n())

Traits<-merge(Traits, tidy_models_TG, by = c("Species"))

#Voltinism (A vs B vs C)
Voltinism <- Traits[Traits$Trait.Group %in% c("A", "B", "C"), ]
kruskal.test(estimate~Trait.Group, data = Voltinism)

#Lecty (A vs D)
Lecty <- Traits[Traits$Trait.Group %in% c("A", "D"), ]
kruskal.test(estimate~Trait.Group, data = Lecty)

#Overwintering Stage (A vs L)
Overwinters <- Traits[Traits$Trait.Group %in% c("A", "L"), ]
kruskal.test(estimate~Trait.Group, data = Overwinters)

#Genus
Traits$Genus <- word(Traits$Species, 1)

Traits %>%
  group_by(Genus) %>%
  summarise(count=n())


Genus <- subset(Traits, Genus == "Andrena" | Genus == "Bombus" | Genus == "Hylaeus" |
                   Genus == "Lasioglossum" | Genus == "Megachile" | Genus == "Osmia" |
                   Genus == "Sphecodes")


kruskal.test(estimate~Genus, data = Genus)
dunnTest(estimate~Genus, data=Genus, method="bh")


#Plot Fig 3
A_cols <- c("#F8766D","#A3A500","#00BF7D")
Voltinism$Voltinism <- paste(Voltinism$Voltinism,"\n\n", sep = "")
Voltinism$Voltinism <- factor(Voltinism$Voltinism , levels=c("Univoltine\n\n", "Bivoltine\n\n", "Variable\n\n"))
A<-ggplot(Voltinism, aes(x = Voltinism, y = estimate, color = Trait.Group))+
  geom_boxplot(outlier.shape = NA)+
  geom_point(position = "jitter")+
  scale_color_manual(values = A_cols) +
  theme_bw()+ 
  theme(axis.text = element_text(size = 10))+
  theme(axis.title = element_text(size = 10, face = "bold"))+
  theme(legend.position = "none")+
  labs(x = "Voltinism", y = "Change in emergence date (Days/°C)", title = expression(bold("Shared Traits")), subtitle= "Period: Spring\nLecty: Polylectic\nOverwintering Stage: Adult within cocoon")+
  ylim(-12, 2)

B_cols <- c("#F8766D", "#00B0F6")
Lecty$Lecty <- paste(Lecty$Lecty,"\n\n", sep = "")
Lecty$Lecty <- factor(Lecty$Lecty, levels=c("Polylectic\n\n", "Oligolectic\n\n"))
B<-ggplot(Lecty, aes(x = Lecty, y = estimate, color = Trait.Group))+
  geom_boxplot(outlier.shape = NA)+
  geom_point(position = "jitter")+
  scale_color_manual(values = B_cols) +
  theme_bw()+ 
  theme(axis.title = element_text(size = 10, face = "bold"))+
  theme(axis.text = element_text(size = 10))+
  theme(legend.position = "none")+
  labs(x = "Lecty", y = "", title = expression(bold("Shared Traits")), subtitle= "Period: Spring\nOverwintering Stage: Adult within cocoon\nVoltinism: Univoltine")+
  ylim(-12, 2)

C_cols <- c("#F8766D", "#E76BF3")
Overwinters$Overwintering.Stage <- factor(Overwinters$Overwintering.Stage, levels=c("Adult within cocoon", "Adult (female only)"))
C<-ggplot(Overwinters, aes(x = Overwintering.Stage, y = estimate, color = Trait.Group))+
  geom_boxplot(outlier.shape = NA)+
  geom_point(position = "jitter")+
  scale_color_manual(values = C_cols) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 6))+
  theme_bw()+
  theme(axis.title = element_text(size = 10, face = "bold"))+
  theme(axis.text = element_text(size = 10))+
  theme(legend.position = "none")+
  labs(x = "Overwintering Stage", y = "", title = expression(bold("Shared Traits")), subtitle= "Period: Spring\nLecty: Polylectic\nVoltinism: Univoltine")+
  ylim(-12, 2)


Legend <- Traits[Traits$Trait.Group %in% c("A", "B", "C", "D", "L"), ]
LegendPlot <- ggplot(Legend, aes(x = Trait.Group, y = estimate, color = Trait.Group))+
  geom_boxplot()+
  geom_point()+
  guides(color = guide_legend(title = "Trait Group"))
Leg <- get_legend(LegendPlot)
Leg <- as_ggplot(Leg)

Blank <-ggplot()+
  theme_void()

D<-ggplot(Genus, aes(x = Genus, y = estimate))+
  geom_boxplot(outlier.shape = NA)+
  geom_point(position = "jitter")+
  theme_bw()+
  labs(x = "Genus", y = "Change in emergence date (Days/°C)")+
  theme(axis.title = element_text(size = 10, face = "bold"))+
  theme(axis.text.x = element_text(face = "italic", size = 10))+
  theme(axis.text.y = element_text(size = 10))+
  theme(legend.position = "none")


Row1 <- ggarrange(A,B,C,Leg, labels = c("A","B","C",""), nrow = 1, widths = c(1,1,1,0.3))
Row2 <- ggarrange(Blank, ncol = 1, labels = c(""))
Row3 <- ggarrange(D, ncol = 1, labels = c("D"))

plot_grid(Row1, Row2, Row3, ncol = 1, rel_heights = c(1, 0.1, 1))