#Analyzing the diet of Alaskan Red knot chicks #Script by Roos Winters #####Start script & load data ##### remove(list=ls()) # clear everything in memory #Load required libraries library(tidyverse) library(cowplot) library(lme4) library(qdapTools) library(writexl) library(readxl) library(devtools) library(vegan) library(mgcv) library(glmmTMB) library(emmeans) library(purrr) #Set working directory setwd("~/Alaskan chick diet - paper/Final_figures/" ) #Make your own colour palette, which ensures that each (taxonomic) group has the same colour in all graphs. custom_palette <- c("Diptera-Anthomyiidae" = "springgreen3", "Diptera-Chironomidae" = "darkseagreen2", "Diptera-Culicidae" = "green3", "Diptera-Empididae" = "darkolivegreen2", "Diptera-Limoniidae" = "seagreen3", "Diptera-Muscidae" = "palegreen4", "Diptera-Phoridae" = "olivedrab3", "Diptera-Sciaridae" = "chartreuse2", "Diptera-Syrphidae" = "green4", "Diptera-Tipulidae" = "darkolivegreen", "Coleoptera-Carabidae" = "#fa9fb5", "Coleoptera-Curculionidae" = "hotpink2", "Coleoptera-Elateridae" = "#dd3497", "Hymenoptera-Braconidae" = "#fdd49e", "Hymenoptera-Ichneumonidae" = "#fc8d59", "Hymenoptera-Tenthredinidae" = "#990000", "Stylommatophora-Agriolimacidae" = "khaki1", "Trombidiformes-Eupodidae" = "purple3", "Lepidoptera-Douglasiidae" = "#efffff", "Lepidoptera-Erebidae" = "#B0E2FA", "Lepidoptera-Geometridae" = "#9ecae1", "Lepidoptera-Noctuidae" = "#3690c0", "Lepidoptera-Nymphalidae" = "#0570b0", "Lepidoptera-Tortricidae" = "#08519c", "Other" = "grey", "2019" = "#117733", "2022" = "#332288", "2023" = "#CC6677", "2024" = "#DDCC77") #we have two datasets: one with the metadata and one with the sequence data --> create one dataset. #Sequence data seq_1 <- read_xlsx("./filtered_family_data_long_new.xlsx") %>% rename(Seq_sample_name = Sample) #Calculate the total RRA (%) per sample --> shows empty samples seq_2 <- seq_1 %>% group_by(Seq_sample_name) %>% summarise(total_diet_perc = sum(RRA)) #Combine the datasets, remove empty samples seq_long <- left_join(seq_1, seq_2, by="Seq_sample_name") %>% filter(!total_diet_perc == 0) #metadata metadata <- read_xlsx("./Mapping_file_all.xlsx")%>% select(!Description) str(metadata$Hatchdate) metadata$Hatchdate <- as.POSIXct(metadata$Hatchdate, origin = "1970-01-01", tz = "UTC") metadata$Hatchdate <- as.Date(metadata$Hatchdate) # convert to Date metadata$collection_date[metadata$collection_date == "N/A"] <- NA metadata$collection_date <- as.Date(as.numeric(metadata$collection_date), origin = "1899-12-30") # Make Year character metadata$Year <- as.character(metadata$Year) #Combine metadata with sequence data md <- left_join(seq_long, metadata, by="Seq_sample_name") %>% arrange(Year, collection_date) %>% mutate(Seq_sample_name = factor(Seq_sample_name, levels = unique(Seq_sample_name))) %>% filter(!Year == "N/A") #The many-to-many relationship is to be expected, as each sample contains multiple taxonomic groups. #Make sure dates are recognized as dates md$collection_date <- as.Date(md$collection_date,format = "%d/%m/%Y") #Change the order of the families (same order are now together) md$taxonomy <- factor(md$taxonomy, levels=c("Coleoptera-Carabidae","Coleoptera-Curculionidae","Coleoptera-Elateridae","Diptera-Anthomyiidae","Diptera-Chironomidae","Diptera-Culicidae","Diptera-Empididae","Diptera-Limoniidae","Diptera-Muscidae","Diptera-Phoridae","Diptera-Sciaridae","Diptera-Syrphidae","Diptera-Tipulidae","Lepidoptera-Douglasiidae","Lepidoptera-Erebidae","Lepidoptera-Geometridae","Lepidoptera-Noctuidae","Lepidoptera-Nymphalidae","Lepidoptera-Tortricidae","Hymenoptera-Braconidae","Hymenoptera-Ichneumonidae","Hymenoptera-Tenthredinidae","Stylommatophora-Agriolimacidae","Trombidiformes-Eupodidae","Other")) snowmelt <- read_xlsx("./snowmelt.xlsx") snowmelt$Year <- as.character(snowmelt$Year) md <- left_join(md, snowmelt, by="Year") #dataframe with only chick data md_chick <- md %>% filter(Age=="Chick") %>% filter(!is.na(taxonomy)) %>% filter(!collection_date == "2022-11-30") %>% filter(!Seq_sample_name %in% c("NIOZ398.013.014","NIOZ400.019.024","NIOZ400.069.074")) %>% mutate(collection_DoY = as.integer(format(collection_date, "%j")), collection_DoY = ifelse(year(collection_date) == 2024, collection_DoY - 1, collection_DoY)) %>% mutate(Hatch_DoY = as.integer(format(Hatchdate, "%j")), Hatch_DoY = ifelse(year(Hatchdate) == 2024, Hatch_DoY - 1, Hatch_DoY)) %>% select(taxonomy, Seq_sample_name, count, RRA, BroodID, RingCode, Year, collection_date, collection_DoY, Hatchdate, Hatch_DoY, Snowmelt_DoY, PlotID) unique(md_chick$Seq_sample_name) md_chick <- md_chick %>% mutate(period = case_when( #2019 collection_date >= as.Date("2019-06-11") & collection_date <= as.Date("2019-06-22") ~ "Early", collection_date >= as.Date("2019-06-23") & collection_date <= as.Date("2019-07-03") ~ "Mid", collection_date >= as.Date("2019-07-04") & collection_date <= as.Date("2019-07-14") ~ "Late", #2022 collection_date >= as.Date("2022-06-16") & collection_date <= as.Date("2022-06-25") ~ "Early", collection_date >= as.Date("2022-06-26") & collection_date <= as.Date("2022-07-04") ~ "Mid", collection_date >= as.Date("2022-07-05") & collection_date <= as.Date("2022-07-13") ~ "Late", # 2023 collection_date >= as.Date("2023-06-20") & collection_date <= as.Date("2023-06-27") ~ "Early", collection_date >= as.Date("2023-06-28") & collection_date <= as.Date("2023-07-05") ~ "Mid", collection_date >= as.Date("2023-07-06") & collection_date <= as.Date("2023-07-12") ~ "Late", # 2024 collection_date >= as.Date("2024-06-16") & collection_date <= as.Date("2024-06-25") ~ "Early", collection_date >= as.Date("2024-06-26") & collection_date <= as.Date("2024-07-04") ~ "Mid", collection_date >= as.Date("2024-07-05") & collection_date <= as.Date("2024-07-13") ~ "Late", TRUE ~ NA_character_)) md_chick$period <- factor(md_chick$period, levels=c("Early", "Mid", "Late")) #Remove all datasets you don't need anymore rm(seq_1, seq_2, metadata, seq_long) #####Chapter 1. Snowmelt and hatchdate##### Snow_hatch <- md_chick %>% select(BroodID, RingCode, Year, Hatchdate, Hatch_DoY, Snowmelt_DoY) %>% unique() %>% group_by(Year, BroodID) %>% mutate(n_chicks = n_distinct(RingCode)) p1 <- Snow_hatch %>% select(RingCode, Year, Hatch_DoY, Snowmelt_DoY) %>% unique() %>% group_by(Year, Hatch_DoY, Snowmelt_DoY) %>% summarise(n_chicks = n_distinct(RingCode)) %>% ggplot(aes(x = Hatch_DoY, y=n_chicks)) + geom_col(aes(fill=Year), color = "black") + geom_vline(aes(xintercept = Snowmelt_DoY), color = "skyblue", linetype = "dashed", linewidth = 1.2) + facet_wrap(~ Year, nrow=4) + scale_y_continuous(breaks = c(0,2,4,6,8,10,12)) + scale_x_continuous(breaks = c(143, 150, 157, 164, 171, 178, 185), labels = c('23 May', "30 May", '06 June', "13 June", "20 June", "27 June", "4 July")) + scale_fill_manual(values = custom_palette) + theme_cowplot() + theme(plot.title = element_text(hjust = 0.5), plot.background = element_rect(fill = "white", color = NA), text = element_text(size = 18), legend.position='none', legend.title=element_blank(), strip.background = element_blank(), strip.placement = "outside", strip.text.y.left = element_text(size = 16, face = "bold"), legend.text = element_text(size = 14), axis.text.x = element_text(size = 16, face = "bold"), axis.title.x = element_blank()) p1 #ggsave("SnowHatch.pdf", plot = phatch, width = 17, height = 8.5, dpi = 300) rm(Snow_hatch, snowmelt) #####Chapter 2. FOO vs. RRA#### FOO_data <- md_chick %>% mutate(presence = ifelse(RRA > 0, 1, 0)) %>% group_by(taxonomy, Seq_sample_name) %>% summarise(presence = max(presence), .groups = "drop") %>% group_by(taxonomy) %>% summarise(Occurrence = sum(presence), .groups = "drop") %>% filter(!taxonomy == "Other") %>% mutate(FOO = Occurrence / 428 * 100) %>% arrange(desc(Occurrence)) %>% mutate(rank=row_number()) pFOO <- FOO_data %>% ggplot(aes(x = taxonomy, y = FOO, fill = taxonomy)) + geom_col(colour = "black", linewidth = 0.5) + labs(y = "FOO (%)") + scale_fill_manual(values = custom_palette) + scale_x_discrete(labels = c("Carabidae", "Curculionidae", "Elateridae", "Anthomyiidae", "Chironomidae", "Culicidae", "Empididae", "Limoniidae", "Muscidae", "Phoridae", "Sciaridae", "Syrphidae","Tipulidae", "Douglasiidae", "Erebidae", "Geometridae", "Noctuidae", "Nymphalidae", "Tortricidae", "Braconidae", "Ichneumonidae","Tenthredinidae", "Agriolimacidae", "Eupodidae")) + #scale_y_continuous(expand = expansion(mult = c(0,0))) + ylim(0,70) + theme_cowplot() + theme(plot.title = element_text(hjust = 0.5), plot.background = element_rect(fill = "white", color = NA), text = element_text(size = 18), legend.position = "none", axis.text.x = element_blank(), axis.title.x = element_blank(), strip.background = element_blank(), strip.text = element_text(size = 16, face = "bold")) RRA_data <- md_chick %>% group_by(taxonomy) %>% summarise(Occurrence = sum(count)) %>% drop_na(c(Occurrence)) %>% filter(!taxonomy == "Other") %>% arrange(desc(Occurrence)) %>% mutate(rank=row_number(), RRA = Occurrence / sum(Occurrence) * 100) RRA_data_samples <- md_chick %>% group_by(Seq_sample_name, taxonomy) %>% summarise(Occurrence = sum(count), .groups = "drop") %>% filter(taxonomy != "Other") %>% group_by(Seq_sample_name) %>% mutate(RRA_sample = Occurrence / sum(Occurrence) * 100) %>% ungroup() %>% group_by(taxonomy) %>% summarise(Avg_RRA = mean(RRA_sample, na.rm = TRUE)) md$taxonomy <- factor(md$taxonomy, levels=c("Coleoptera-Carabidae","Coleoptera-Curculionidae","Coleoptera-Elateridae","Diptera-Anthomyiidae","Diptera-Chironomidae","Diptera-Culicidae","Diptera-Empididae","Diptera-Limoniidae","Diptera-Muscidae","Diptera-Phoridae","Diptera-Sciaridae","Diptera-Syrphidae","Diptera-Tipulidae","Lepidoptera-Douglasiidae","Lepidoptera-Erebidae","Lepidoptera-Geometridae","Lepidoptera-Noctuidae","Lepidoptera-Nymphalidae","Lepidoptera-Tortricidae","Hymenoptera-Braconidae","Hymenoptera-Ichneumonidae","Hymenoptera-Tenthredinidae","Stylommatophora-Agriolimacidae","Trombidiformes-Eupodidae","Other")) pRRA <- RRA_data_samples %>% mutate(taxonomy = factor(taxonomy, levels = taxonomy)) %>% ggplot(aes(x = taxonomy, y = Avg_RRA, fill = taxonomy)) + geom_col(colour = "black", linewidth = 0.5) + labs(y = "RRA (%)", x = NULL) + ylim(0, 50) + scale_fill_manual(values = custom_palette) + #scale_y_continuous(expand = expansion(mult = c(0,0))) + scale_x_discrete(labels = c("Carabidae", "Curculionidae", "Elateridae", "Anthomyiidae", "Chironomidae", "Culicidae", "Empididae", "Limoniidae", "Muscidae", "Sciaridae", "Syrphidae","Tipulidae", "Douglasiidae", "Erebidae", "Geometridae", "Noctuidae", "Nymphalidae", "Tortricidae", "Braconidae", "Ichneumonidae","Tenthredinidae", "Agriolimacidae", "Eupodidae")) + theme_cowplot() + theme( plot.title = element_text(hjust = 0.5), plot.background = element_rect(fill = "white", color = NA), text = element_text(size = 18), legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), axis.title.x = element_blank(), strip.background = element_blank(), strip.text = element_text(size = 16, face = "bold") ) p2 <- cowplot::plot_grid( pFOO, pRRA, labels = c("A", "B"), ncol = 1, nrow = 2, align = "v",axis = "tb",rel_heights = c(1, 1.2)) p2 #ggsave("FOOvsRRA.pdf", plot = p2, width = 14, height = 8.5, dpi = 300) #####Chapter 3: Diet comparison between years ##### #calculate the yearly means & average of the means of 2022 and 2023 year_comparison <- md_chick %>% group_by(Year, taxonomy) %>% summarise(Occurrence_year = sum(RRA)) %>% drop_na(c(Year, Occurrence_year)) %>% ungroup() %>% group_by(Year) %>% mutate(Total_reads_year = sum(Occurrence_year)) %>% mutate(percentage_year = (Occurrence_year/Total_reads_year)*100) %>% ungroup() %>% group_by(taxonomy) %>% mutate(mean_percentage = mean(percentage_year)) #number of samples per year nr_samples <- md_chick %>% select(c(Year, Seq_sample_name)) %>% unique() %>% group_by(Year) %>% summarise(n = n()) %>% filter(!Year== "N/A") %>% drop_na() #plot graph comparing the chick diet of 2019, 2022, 2023 and 2024 p3 <- year_comparison %>% ggplot(aes(x=Year, y=percentage_year, fill=taxonomy)) + geom_bar(stat="identity", colour ="black", linewidth = 0.5) + labs(y = "RRA (%)") + scale_fill_manual(values = custom_palette, guide = guide_legend(ncol = 1))+ scale_fill_manual(values = custom_palette, guide = guide_legend(ncol = 1), labels = c("Carabidae", "Curculionidae", "Elateridae", "Anthomyiidae", "Chironomidae", "Culicidae", "Empididae", "Limoniidae", "Muscidae", "Phoridae", "Sciaridae", "Syrphidae","Tipulidae", "Douglasiidae", "Erebidae", "Geometridae", "Noctuidae", "Nymphalidae", "Tortricidae", "Hydrobiosidae", "Braconidae", "Ichneumonidae", "Pteromalidae", "Tenthredinidae", "Littorinidae", "Agriolimacidae", "Succineidae", "Ixodidae", "Eupodidae", "Tetranychidae", "Other")) + geom_text(data = nr_samples, aes(x = Year, y = 105, label = paste0(n)),inherit.aes = FALSE, size = 4, vjust = 1) + scale_y_continuous(expand = expansion(mult = c(0,0))) + theme_cowplot() + theme(plot.title = element_text(hjust = 0.5), plot.background = element_rect(fill = "white", color = NA), text = element_text(size = 18), legend.position='right', legend.title=element_blank(), axis.title.x = element_blank(), axis.text.x = element_text(vjust = 0.5), strip.background = element_blank(), strip.text = element_text(size = 16, face = "bold"), legend.text = element_text(size = 14)) p3 #save plot & save dataset #ggsave("yearly_average_diet.pdf", plot = p3, width = 8.5, height = 12, dpi = 300) #write_xlsx(md,"./combined_data.xlsx") #Mean and SE for all samples from 2019 summary_2019 <- md_chick %>% filter(Year == 2019) %>% group_by(taxonomy) %>% summarise(mean_percentage_2019 = mean(RRA, na.rm = TRUE), st_dev_2019 = sd(RRA, na.rm = TRUE), st_error_2019 = st_dev_2019 / sqrt(n()), .groups = "drop") #Mean and SE for all samples from 2022 summary_2022 <- md_chick %>% filter(Year == 2022) %>% group_by(taxonomy) %>% summarise(mean_percentage_2022 = mean(RRA, na.rm = TRUE), st_dev_2022 = sd(RRA, na.rm = TRUE), st_error_2022 = st_dev_2022 / sqrt(n()), .groups = "drop") #Mean and SE for all samples from 2023 summary_2023 <- md_chick %>% filter(Year == 2023) %>% group_by(taxonomy) %>% summarise(mean_percentage_2023 = mean(RRA, na.rm = TRUE), st_dev_2023 = sd(RRA, na.rm = TRUE), st_error_2023 = st_dev_2023 / sqrt(n()), .groups = "drop") #Mean and SE for all samples from 2024 summary_2024 <- md_chick %>% filter(Year == 2024) %>% group_by(taxonomy) %>% summarise(mean_percentage_2024 = mean(RRA, na.rm = TRUE), st_dev_2024 = sd(RRA, na.rm = TRUE), st_error_2024 = st_dev_2024 / sqrt(n()), .groups = "drop") #Merge summary_combined <- summary_2022 %>% left_join(summary_2023, by = "taxonomy") %>% left_join(summary_2019, by = "taxonomy") %>% left_join(summary_2024, by = "taxonomy") %>% mutate( mean_percentage_combined = (mean_percentage_2019 + mean_percentage_2022 + mean_percentage_2023 + mean_percentage_2024) / 4, st_error_combined = sqrt((st_error_2019^2 + st_error_2022^2 + st_error_2023^2 + st_error_2024^2) / 4) ) %>% mutate(across(where(is.numeric), ~ round(.x, 1))) %>% select(taxonomy,mean_percentage_2019, st_error_2019, mean_percentage_2022, st_error_2022, mean_percentage_2023, st_error_2023, mean_percentage_2024, st_error_2024, mean_percentage_combined) #write_xlsx(summary_combined,"./occurrance_table.xlsx" ) rm(nr_samples, summary_2022, summary_2023, summary_combined, year_comparison) #####Chapter 4: Seasonal variation #### #change in diet over the season (early, mid and late season) nr_samples <- md_chick %>% select(c(Year, collection_DoY, Seq_sample_name)) %>% unique() %>% group_by(Year, collection_DoY) %>% summarise(n = n()) %>% drop_na() %>% filter(!Year == "N/A") #Calculate the average diet per day daily_diet <- md_chick %>% group_by(Year, collection_DoY, taxonomy) %>% summarise(Occurrence_day = sum(RRA)) %>% ungroup() %>% group_by(Year, collection_DoY) %>% mutate(Total_reads_day = sum(Occurrence_day)) %>% mutate(percentage_day = (Occurrence_day/Total_reads_day)*100) %>% ungroup() #Graph one bar per day p4 <- daily_diet %>% ggplot(aes(x = collection_DoY, y = percentage_day, fill = taxonomy)) + geom_bar(stat = "identity", colour = "black", linewidth = 0.5) + labs(y = "% Diet", x = "Date") + scale_fill_manual(values = custom_palette, guide = guide_legend(ncol = 1), labels = c("Carabidae", "Curculionidae", "Elateridae", "Anthomyiidae", "Chironomidae", "Culicidae", "Empididae", "Limoniidae", "Muscidae", "Phoridae", "Sciaridae", "Syrphidae", "Tipulidae","Douglasiidae", "Erebidae", "Geometridae", "Noctuidae", "Nymphalidae", "Tortricidae", "Hydrobiosidae", "Braconidae", "Ichneumonidae","Pteromalidae", "Tenthredinidae", "Littorinidae", "Agriolimacidae", "Succineidae", "Ixodidae", "Eupodidae", "Tetranychidae", "Other")) + geom_text(data = nr_samples, aes(x = collection_DoY, y = 100, label = n), inherit.aes = FALSE, size = 4,vjust = 0,nudge_y = 2) + scale_x_continuous(breaks = c(170, 175, 180, 185, 190, 195), labels = c('19 June', '24 June', '29 June', '04 July', '09 July', '14 July')) + facet_grid(vars(Year)) + theme_cowplot() + theme(plot.title = element_text(hjust = 0.5), text = element_text(size = 18), legend.position = "right", legend.title = element_blank(), strip.background = element_blank(), strip.text = element_text(size = 16, face = "bold"), legend.text = element_text(size = 14), axis.title.x = element_blank()) + coord_cartesian(ylim = c(0, 100), clip = "off") p4 #ggsave("daily_diet.png", plot = p4, width = 17, height = 14, dpi = 300) #write_xlsx(daily_diet, "./Daily_diet_chicks.xlsx") #Divide the data into three (~ equal) periods periods <- md_chick %>% group_by(Year, period, taxonomy) %>% summarise(Occurrence_period = sum(RRA)) %>% drop_na() %>% ungroup() %>% group_by(Year, period) %>% mutate(Total_reads_period = sum(Occurrence_period)) %>% mutate(percentage_period = (Occurrence_period/Total_reads_period)*100) %>% ungroup() nr_samples_periods <- md_chick %>% distinct(Year, period, Seq_sample_name) %>% # force unique samples group_by(Year, period) %>% summarise(n = n(), .groups = "drop") %>% drop_na() p5 <- periods %>% ggplot(aes(x=period, y=percentage_period, fill=taxonomy)) + geom_bar(stat="identity", colour ="black", linewidth = 0.5) + geom_text(data = nr_samples_periods, aes(x = period, y = 110, label = paste0(n)),inherit.aes = FALSE, size = 4, vjust = 1) + labs(y = "RRA (%)", x= "Period") + scale_y_continuous(breaks = seq(0, 100, by = 25)) + scale_x_discrete(breaks = c("Early", "Mid", "Late")) + scale_fill_manual(values = custom_palette, guide = guide_legend(ncol = 1), labels = c("Carabidae","Curculionidae","Elateridae", "Anthomyiidae", "Chironomidae", "Culicidae", "Empididae","Limoniidae","Muscidae", "Sciaridae","Syrphidae","Tipulidae", "Douglasiidae","Erebidae","Geometridae","Noctuidae", "Nymphalidae", "Tortricidae","Braconidae","Ichneumonidae", "Tenthredinidae", "Agriolimacidae","Eupodidae","Other")) + facet_wrap(vars(Year), ncol=1, strip.position = "left") + theme_cowplot() + theme(plot.title = element_text(hjust = 0.5), plot.background = element_rect(fill = "white", color = NA), text = element_text(size = 18), legend.position='right', legend.title=element_blank(), strip.background = element_blank(), strip.placement = "outside", strip.text.y.left = element_text(size = 16, face = "bold"),, legend.text = element_text(size = 14), axis.text.x = element_text(size = 16, face = "bold"), axis.title.x = element_blank()) p5 #ggsave("periods_new.pdf", plot = p5, width = 8.5, height = 12, dpi = 300) #write_xlsx(md, "./combined_data.xlsx") rm(periods, nr_samples, nr_samples_periods) #FOO per period FOO_data <- md_chick %>% mutate(presence = ifelse(RRA > 0, 1, 0)) %>% group_by(taxonomy, Seq_sample_name, period, Year) %>% summarise(presence = max(presence), .groups = "drop") %>% group_by(taxonomy, Year, period) %>% summarise(Occurrence = sum(presence), .groups = "drop") %>% filter(taxonomy != "Other") %>% left_join(md_chick %>% distinct(Seq_sample_name, period, Year) %>% count(period, Year, name = "n_samples"), by = c("period", "Year")) %>% mutate(FOO = Occurrence / n_samples * 100) %>% arrange(Year, period, desc(Occurrence)) %>% group_by(Year, period) %>% mutate(rank = row_number()) %>% ungroup() p6 <- ggplot(FOO_data,aes(x = forcats::fct_rev(taxonomy), y = FOO, fill = taxonomy)) + geom_col(colour = "black", linewidth = 0.4) + facet_grid(Year ~ period, scales = "free_x",space = "free_x") + scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 25), expand = expansion(mult = c(0, 0.05))) + scale_fill_manual(values = custom_palette, guide = guide_legend(ncol = 1), labels = c("Carabidae","Curculionidae","Elateridae", "Anthomyiidae", "Chironomidae", "Culicidae", "Empididae","Limoniidae","Muscidae", "Sciaridae","Syrphidae","Tipulidae", "Douglasiidae","Erebidae","Geometridae","Noctuidae", "Nymphalidae", "Tortricidae","Braconidae","Ichneumonidae", "Tenthredinidae", "Agriolimacidae","Eupodidae","Other")) + labs(y = "FOO (%)", x = NULL) + theme_cowplot() + theme(text = element_text(size = 16), axis.text.x.top = element_blank(), axis.ticks.x = element_blank(), axis.text.y = element_text(size = 14), strip.text = element_text(size = 14, face = "bold"), strip.background = element_blank(), legend.position = "none") p6 #ggsave("periods_FOO.pdf", plot = p6, width = 10, height = 8.5, dpi = 300) md_chick <- md_chick %>% mutate(Age_chick = as.numeric(collection_date - Hatchdate)) taxa_list <- c("Coleoptera-Curculionidae","Diptera-Tipulidae", "Lepidoptera-Noctuidae","Hymenoptera-Ichneumonidae", "Stylommatophora-Agriolimacidae") #DoY models fit_taxon_model_DoY <- function(taxon_name, df) { taxon_data <- df %>% filter(taxonomy == taxon_name, Year != 2019) %>% mutate(RRA_prop_raw = RRA / 100, n_obs = n(), RRA_prop = (RRA_prop_raw * (n_obs - 1) + 0.5) / n_obs, Year = factor(Year)) model <- glmmTMB(RRA_prop ~ collection_DoY * Year, family = beta_family(link = "logit"),data = taxon_data) slopes <- emtrends(model, specs = "Year", var = "collection_DoY") %>% test() %>% as.data.frame() %>% mutate(Taxon = taxon_name, Temporal_scale = "DoY") list(model = model, slopes = slopes, data = taxon_data) } results_DoY <- map(taxa_list, fit_taxon_model_DoY, df = md_chick) slopes_DoY <- map_dfr(results_DoY, "slopes") slopes_DoY #Age models fit_taxon_model_Age <- function(taxon_name, df) { taxon_data <- df %>% filter(taxonomy == taxon_name, Year != 2019) %>% mutate(RRA_prop_raw = RRA / 100, n_obs = n(), RRA_prop = (RRA_prop_raw * (n_obs - 1) + 0.5) / n_obs, Year = factor(Year)) model <- glmmTMB( RRA_prop ~ Age_chick * Year, family = beta_family(link = "logit"), data = taxon_data) slopes <- emtrends(model, specs = "Year", var = "Age_chick") %>% test() %>% as.data.frame() %>% mutate(Taxon = taxon_name, Temporal_scale = "Age") list(model = model, slopes = slopes, data = taxon_data) } results_Age <- map(taxa_list, fit_taxon_model_Age, df = md_chick) slopes_Age <- map_dfr(results_Age, "slopes") #combine them slopes_all <- left_join(slopes_DoY, slopes_Age, by=c("Year", "Taxon")) #Save as CSV for easy download #write.csv(slopes_all, "slopes_DoY_Age.csv", row.names = FALSE) #####Chapter 5: Brood variation ##### md_37M_2023 <- md_chick %>% filter(Year == "2023") %>% filter(PlotID == "37M") %>% mutate(BroodID = ifelse(BroodID == "1MT/1CX", "1MT", BroodID)) %>% filter(BroodID %in% c("1MP", "1KU", "1MH", "1CX", "1AY", "1MT","1NL")) md_37M_2024 <- md_chick %>% filter(Year == "2024") %>% filter(PlotID == "37M") %>% filter(BroodID %in% c("1MP", "1KU", "1MT","1NL")) full_dates <- seq(as.Date("2023-06-19"), as.Date("2023-07-11"), by = "day") ranks_2023 <- md_37M_2023 %>% select(collection_date, BroodID, RingCode, Seq_sample_name) %>% unique() %>% arrange(BroodID, collection_date, RingCode) %>% group_by(collection_date, BroodID) %>% mutate(rank=row_number()) ranks_2024 <- md_37M_2024 %>% select(collection_date, BroodID, RingCode, Seq_sample_name) %>% unique() %>% arrange(BroodID, collection_date, RingCode) %>% group_by(collection_date, BroodID) %>% mutate(rank=row_number()) data_37M_2023 <- md_37M_2023 %>% left_join(ranks_2023, by = c("collection_date", "BroodID", "RingCode", "Seq_sample_name")) %>% filter(!Seq_sample_name %in% c("NIOZ398.163.164", "NIOZ399.145.148", "NIOZ399.013.016")) data_37M_2024 <- md_37M_2024 %>% left_join(ranks_2024, by = c("collection_date", "BroodID", "RingCode", "Seq_sample_name")) duplos_2023 <- data_37M_2023 %>% select(c(BroodID, collection_date, RingCode, Seq_sample_name, rank)) %>% group_by(BroodID, collection_date, RingCode) %>% unique() %>% mutate(n = n()) %>% filter(n >= 2) %>% mutate(duplicate = ifelse(collection_date == "2023-07-05" & RingCode == "1232-33635", "^", "*")) duplos_2024 <- data_37M_2024 %>% select(c(BroodID, collection_date, RingCode, Seq_sample_name, rank)) %>% group_by(BroodID, collection_date, RingCode) %>% unique() %>% mutate(n = n()) %>% filter(n >= 2) %>% mutate(duplicate = "*") data_37M_2023 <- left_join(data_37M_2023, duplos_2023, by=c("BroodID", "collection_date", "RingCode", "Seq_sample_name", "rank")) data_37M_2024 <- left_join(data_37M_2024, duplos_2024, by=c("BroodID", "collection_date", "RingCode", "Seq_sample_name", "rank")) p2023 <- data_37M_2023 %>% ggplot(aes(x = rank, y = RRA, fill = taxonomy)) + geom_bar(stat = "identity", position = "stack", col = "black", linewidth = 0.5, width = 0.9) + geom_text(aes(x = rank, y = -6, label = duplicate), inherit.aes = FALSE, size = 3, color = "black") + labs(y = "% Diet") + scale_fill_manual(values = custom_palette, guide = guide_legend(ncol = 1)) + scale_y_continuous(expand = expansion(mult = c(0, 0))) + facet_grid(factor(BroodID, levels=c('1MP','1NL','1MH', '1CX', "1MT", "1KU")) ~ factor(collection_date), scales = "free_x", space = "free", labeller = as_labeller(c("2023-06-24" = "24 June", "2023-06-26" = "26 June", "2023-06-29" = "29 June", "2023-07-02" = "02 July", "2023-07-05" = "05 July", "2023-07-08" = "08 July", "2023-07-11" = "11 July", "1CX" = "1CX", "1KU" = "1KU", "1MH" = "1MH", "1MP" = "1MP", "1MT" = "1MT", "1NL" = "1NL"))) + theme_cowplot() + theme(plot.title = element_text(hjust = 0.5), plot.background = element_rect(fill = "white", color = NA), text = element_text(size = 18), legend.position = 'none', legend.title=element_blank(), axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(), strip.background = element_blank(), strip.text.y = element_text(size = 16, face = "bold", angle = 0), panel.spacing.x = unit(0.5, "lines")) p2023 p2024 <- data_37M_2024 %>% ggplot(aes(x = rank, y = RRA, fill = taxonomy)) + geom_bar(stat = "identity", position = "stack", col = "black", linewidth = 0.5, width = 0.9) + geom_text(aes(x = rank, y = -6, label = duplicate), inherit.aes = FALSE, size = 3, color = "black") + labs(y = "% Diet") + scale_fill_manual(values = custom_palette) + scale_y_continuous(expand = expansion(mult = c(0, 0))) + facet_grid(factor(BroodID, levels=c('1MT','1KU','1MP', '1NL')) ~ factor(collection_date), space = "fixed", labeller = as_labeller(c("2024-06-22" = "22 June", "2024-06-24" = "24 June", "2024-06-25" = "25 June","2024-06-26" = "26 June", "2024-06-27" = "27 June","2024-06-29" = "29 June","2024-07-01" = "01 July", "2024-07-04" = "04 July","2024-07-06" = "06 July","2024-07-09" = "09 July","2024-07-11" = "11 July", "1KU" = "1KU", "1MP" = "1MP", "1MT" = "1MT", "1NL" = "1NL"))) + theme_cowplot() + theme(plot.title = element_text(hjust = 0.5), plot.background = element_rect(fill = "white", color = NA), text = element_text(size = 18), legend.position = 'none', legend.title=element_blank(), axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(), strip.background = element_blank(), strip.text.y = element_text(size = 16, face = "bold", angle = 0), panel.spacing.x = unit(0.5, "lines")) p2024 make_plot <- function(df, year_label) { ggplot(df, aes(x = factor(rank), y = RRA, fill = taxonomy)) + geom_bar(stat = "identity", position = "stack", col = "black", linewidth = 0.5, width = 0.8) + geom_text(aes(x = factor(rank), y = -6, label = duplicate), inherit.aes = FALSE, size = 3, color = "black") + labs(y = "% Diet", title = year_label) + scale_fill_manual(values = custom_palette, guide = guide_legend(ncol = 1)) + scale_y_continuous(expand = expansion(mult = c(0, 0))) + facet_wrap(factor(BroodID, levels = c('1MP','1NL','1MH','1CX','1MT','1KU')) ~ facet_label,scales = "free_x",space = "fixed") + theme_cowplot() + theme( plot.title = element_text(hjust = 0.5), text = element_text(size = 18), legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(), strip.background = element_blank(), strip.text.y = element_text(size = 16, face = "bold", angle = 0), panel.spacing.x = unit(0.5, "lines")) } p2023 <- make_plot(data_37M_2023, "2023") p2024 <- make_plot(data_37M_2024, "2024") library(patchwork) p7 <- p2023 / p2024 + plot_layout(heights = c(6, 4), guides = "collect") p7 #ggsave("brood_variation_newcols.pdf", plot = p7, width = 17, height = 12, dpi = 300) #####Chapter 6: Pitfall data ##### #Load data and data prep pitfall <- read_xlsx("./pitfalls_2022_23_24.xlsx")%>% mutate(plotID = ifelse(Station == "Grid42" | Station == "Grid45", "37M", plotID)) %>% rename(Plot.ID = plotID) %>% filter(!Poskey %in% c("AK22_076","AK22_077", "AK22_078", "AK22_079", "AK22_080", "AK22_081", "AK22_082", "AK22_083", "AK22_084", "AK22_085", "AK22_086", "AK22_087", "AK22_088", "AK22_089", "AK22_090")) check <- pitfall %>% filter(Year == 2022) %>% filter(Date == 179) pitfall$Date <- as.Date(pitfall$Date ,format = "%m/%d/%y") pitfall$Date <- as.numeric(format(pitfall$Date, "%j")) pitfall$Station[pitfall$Station == "Grid7"] <- "Grid07" pitfall$Plot.ID[pitfall$Station == "Grid07"] <- "37M" pitfall2 <- pitfall %>% select(Year, Date, Plot.ID, Station) %>% unique() %>% group_by(Year, Date) %>% mutate(nr_pitfalls_day = n()) %>% ungroup() all_combinations <- pitfall2 %>% distinct(Year, Date, Plot.ID) %>% expand_grid(taxonomy = unique(pitfall$Family)) %>% filter(taxonomy %in% c("Curculionidae", "Tipulidae", "Elateridae")) pitfall3 <- left_join(all_combinations, pitfall2) %>% mutate(Biomass = 0) #many-to-many relationship is expected pitfall4 <- pitfall %>% select(Year, Date, Plot.ID, Family, Biomass) %>% filter(Family %in% c("Curculionidae", "Tipulidae", "Elateridae")) %>% rename(taxonomy = Family) pitfall_fam <- left_join(pitfall3, pitfall4, by=c("Year", "Date", "Plot.ID", "taxonomy")) %>% mutate(Biomass.x = ifelse(is.na(Biomass.y), Biomass.x, Biomass.y)) %>% select(!Biomass.y) %>% rename(Biomass = Biomass.x) #many-to-many is expected pitfall_fam <- pitfall_fam %>% group_by(Year, Date, taxonomy) %>% mutate(total_biomass = sum(Biomass, na.rm = TRUE), daily_biomass = (total_biomass/5), average_daily_biomass = (daily_biomass/nr_pitfalls_day))%>% ungroup() %>% select(!c(Biomass, Plot.ID, Station)) %>% unique() %>% rename(taxonomy = taxonomy) %>% mutate(taxonomy = case_when(taxonomy == "Curculionidae" ~ "Coleoptera-Curculionidae", taxonomy=="Tipulidae" ~ "Diptera-Tipulidae", taxonomy=="Elateridae" ~ "Coleoptera-Elateridae")) %>% filter(!taxonomy == "Coleoptera-Elateridae") rm(pitfall2,pitfall3,pitfall4, all_combinations) #####Combining pitfall data and family data family_data <- md_chick %>% rename(Date = collection_date) %>% filter(taxonomy == "Coleoptera-Curculionidae" | taxonomy == "Diptera-Tipulidae") %>% filter(!Year == "2019") family_data$Date <- as.Date(family_data$Date ,format = "%m/%d/%y") family_data$Date <- as.numeric(format(family_data$Date , "%j")) family_data$Year <- as.character(family_data$Year) coeff <- 5 chick_periods <- md_chick %>% filter(!Year==2019) %>% group_by(Year) %>% summarise(start_chick = min(Hatch_DoY), end_chick = max(collection_DoY)) %>% filter(Year %in% c(2022, 2023, 2024)) p8 <- ggplot() + geom_rect(data = chick_periods, aes(xmin = start_chick, xmax = end_chick, ymin = -Inf, ymax = Inf), fill = "snow4", alpha = 0.15, inherit.aes = FALSE) + geom_bar(data = pitfall_fam, aes(x = Date, y = average_daily_biomass, fill=taxonomy), col= "black", stat = "identity", size = .1) + #geom_text(data=pitfall_fam, aes(x = Date, y = average_daily_biomass, label = ifelse(average_daily_biomass == 0,"*", "")), vjust = -0.5, color = "black", size=3) + geom_vline(data = family_data, aes(xintercept = Snowmelt_DoY), color = "skyblue", linetype = "dashed", linewidth = 1.2) + scale_y_continuous(name = "Biomass (mg/pitfall/day)", limits = c(0, 20))+ scale_x_continuous(breaks = c(143, 150, 157, 164, 171, 178, 185, 192, 199, 206),labels = c('23 May','30 May','06 June','13 June','20 June','27 June','04 July','11 July', '18 July', '25 July')) + facet_grid(Year ~ .) + scale_fill_manual(values = custom_palette) + scale_colour_manual(values = custom_palette) + theme_cowplot() + theme(plot.title = element_text(hjust = 0.5), text = element_text(size = 18), legend.position = 'none', axis.title.x = element_blank(), strip.text = element_text(size = 16, face = "bold"), strip.placement = "outside", strip.background = element_blank()) p8 #ggsave("Pitfall_diet.pdf", plot = p8, width = 17, height = 8.5, dpi = 300) #write_xlsx(md_chick, "./Rscripts/metadata_chicks.xlsx")