# CHAPTER 1 # 1. INSTALL AND LOAD PACKAGES -------------------------------------------- remotes::install_github("ropensci/rnaturalearthhires") library(broom) library(codyn) library(ggplot2); theme_set(theme_bw()) library(ggspatial) library(iNEXT) library(knitr) library(patchwork) library(RColorBrewer) library(rnaturalearth) library(rnaturalearthdata) library(rnaturalearthhires) library(scatterpie) library(sf) library(terra) library(tidyverse) library(vegan) options(stringsAsFactors = F) # 2. MAP OF THE STUDY AREA ------------------------------------------------ # 2.1. LOAD DATA AND SHAPEFILES -------------------------------------- # Study site data data <- read.csv("Chapter1_data/proporcao.csv", h = T) site <- data %>% # Sum mutate(TOTAL = rowSums(.[5:15])) %>% # Proportion mutate(multiplier = log10(TOTAL) / log10(max(TOTAL))) site # Brazil (medium resolution) br <- rnaturalearth::ne_countries(scale = "medium", country = "brazil", returnclass = "sf") # Brazilian states states <- rnaturalearth::ne_states(country = "brazil", returnclass = "sf") # Brazilian biomes urlRemote <- "http://geoftp.ibge.gov.br/" pathRemote <- "informacoes_ambientais/estudos_ambientais/biomas/vetores/" fileName <- "Biomas_250mil.zip" url <- paste0(urlRemote, pathRemote, fileName) file <- basename(url) download.file(url, file) # Unzip biomes_unzip <- unzip("Biomas_250mil.zip", exdir = "Chapter1_data/Biomes_250mil") biomes_unzip # Load vector data biomes <- terra::vect("Chapter1_data/Biomes_250mil/lm_bioma_250.shp") # Convert to sf biomes_sf <- sf::st_as_sf(biomes) biomes_sf # Cerrado polygon cerrado <- subset(biomes_sf, Bioma == "Cerrado") # Amazonia polygon amazonia <- subset(biomes_sf, Bioma == "Amazônia") # Transition polygons forest <- terra::vect("Chapter1_data/Shapes/Disjunct_forest_matrix/", "disjunt_matrix_forest") savanna <- terra::vect("Chapter1_data/Shapes/Disjunct_savanna_matrix/", "disjunt_matrix_savana") # Convert to sf forest_sf <- sf::st_as_sf(forest) savanna_sf <- sf::st_as_sf(savanna) # 2.2. CREATE MAP ---------------------------------------------------- # Map fig1 <- ggplot() + # Brazil geom_sf(data = br, fill = NA, colour = "dark grey", linewidth = 0.5) + # Brazilian states geom_sf(data = states, fill = NA, colour = "dark grey", linewidth = 0.2, show.legend = F) + # Cerrado geom_sf(data = cerrado, fill = "#d73027", alpha = 0.2, colour = "#d73027", linewidth = 0.1) + # Amazonia geom_sf(data = amazonia, fill = "#1a9850", alpha = 0.2, colour = "#1a9850", linewidth = 0.1) + # Transition: savanna geom_sf(data = savanna_sf, fill = "#d73027", alpha = 0.6, colour = "#d73027", linewidth = 0.2) + # Transition: forest geom_sf(data = forest_sf, fill = "#1a9850", alpha = 0.6, colour = "#1a9850", linewidth = 0.2) + # Study site geom_point(data = site, aes(Longitude, Latitude), shape = 21, size = 3, fill = "white", stroke = 2) + # Axes labels xlab("") + ylab("") + # Scale ggspatial::annotation_scale(location = "bl", width_hint = 0.3) + # North arrow ggspatial::annotation_north_arrow(location = "bl", which_north = "true", height = unit(0.5, "in"), width = unit(0.5, "in"), pad_x = unit(0.2, "in"), pad_y = unit(0.5, "in"), style = north_arrow_fancy_orienteering) + # Coordinate system coord_sf() # View and sava fig1 ggsave("Chapter1_data/Figures/Figure1.pdf", plot = fig1, width = 8, height = 8) # 3. RAREFACTION ANALYSIS ------------------------------------------------- RareCamp <- read.table ("Chapter1_data/rarefy.txt", row.names = 1, h = T) Ext.Camp <- iNEXT::iNEXT(RareCamp, datatype = "abundance", se = T, conf = 0.95, endpoint = NULL, nboot = 100) # Basic info Ext.Camp$DataInfo # Diversity asymptotic estimates Ext.Camp$AsyEst # Richness estimates Ext.Camp$iNextEst # Graphs fig4a <- ggiNEXT(Ext.Camp, type = 1, color.var = "Order.q") + theme_bw(base_size = 16) + theme(legend.position = "bottom") fig4a ggsave("Chapter1_data/Figures/Figure4A.pdf", plot = fig4a, width = 8, height = 8) fig4b <- ggiNEXT(Ext.Camp, type = 2, color.var = "Order.q") + theme_bw(base_size = 16) + theme(legend.position = "bottom") fig4b ggsave("Chapter1_data/Figures/Figure4B.pdf", plot = fig4b, width = 8, height = 8) # 4. TEMPORAL COMMUNITY DYNAMICS ------------------------------------------ # Read data snakec <- read.csv("Chapter1_data/data_snakes.csv", sep= ';', h = T) head(snakec) # Species abundances totais <- aggregate(abundance ~ species, sum, data = snakec) # Ordered abundances totais_ord <- totais[order(totais$abundance, decreasing = T),] totais_ord # Add rank column totais_ord$rank <- seq_along(totais_ord$abundance) # Graph plot(totais.ord$abundance, pch = 19, type = "b", lwd = 0.5, las = 1, bty = "n", ylab = "Abundance", xlab = "Rank") # Plot Rank-Abundance Distribution (RAD) with all species labels rad_plot <- ggplot(totais_ord, aes(x = rank, y = abundance)) + # Bar plot geom_bar(stat = "identity", fill = "skyblue") + # Add species labels geom_text(aes(label = species), # Position labels to the right of the bars hjust = -0.1, # Center vertically along the bars vjust = 0.5, # Adjust label text size size = 3) + labs(x = "Rank", y = "Abundance") + # Flip coordinates for better label placement coord_flip() + theme_minimal() + # Remove legend if not needed theme(legend.position = "none") # Display the plot rad_plot # 4.1. TEMPORAL DIVERSITY PATTERNS ----------------------------------- # 4.1.1. VISUALIZING TEMPORAL DYNAMICS WITH RANK CLOCKS --------- # Define the species of interest species_of_interest <- c("pseudoboa nigra", "erythrolamprus poecilogyrus", "sibynomorphus mikanii") # Create a complete grid of all seasons and species complete_grid <- expand.grid(season = unique(snakec$season), species = species_of_interest) # Aggregate the data aggdat <- aggregate(abundance ~ season + species, data = subset(snakec, species %in% species_of_interest), FUN = sum) # Merge the aggregated data with the complete grid, filling missing values with zero aggdat_complete <- merge(complete_grid, aggdat, by = c("season", "species"), all.x = TRUE) # Replace NA values with 0 in the abundance column aggdat_complete$abundance[is.na(aggdat_complete$abundance)] <- 0 # View the result aggdat_complete # Rename species for graphing purposes aggdat$species <- as.factor(as.character(aggdat$species)) levels(aggdat$species) = c("E. poecilogyrus", "P. nigra", "S. mikanii") # Create the graph rankshift.graph <- ggplot(aggdat, aes(season, abundance, color = species)) + geom_line(linewidth = 0.5) + coord_polar() + labs(x = "Season", y = "Abundance", color = "Species") + scale_x_continuous(breaks = 1:6, labels = c("Dry 2018", "Wet 2018", "Dry 2019", "Wet 2019", "Dry 2020", "Wet 2020")) + scale_color_manual(values = c("#009E73", "#e79f00", "#9ad0f3")) + geom_segment(aes(x = 1, y = 0, xend = 1, yend = 10), color = "grey70") + theme_bw() + theme(text = element_text(size = 14), strip.text.x = element_text(size = 14), strip.background = element_blank(), panel.grid.major = element_line(size = 1)) + theme(legend.position = "bottom", legend.text = element_text(face = "italic")) rankshift.graph # 4.2. QUANTIFYING TEMPORAL DYNAMICS --------------------------------- # 4.2.1. SPECIES RICHNESS --------------------------------------- # Calculate richness within seasons rich.dat <-snakec %>% filter(abundance > 0) %>% mutate(rich = 1) %>% group_by(season) %>% summarize(richness = sum(rich)) %>% tibble::as_tibble() # Create the graph rich.graph <- ggplot(rich.dat, aes(x = season, y = richness)) + geom_line(linewidth = 1) + scale_x_continuous(breaks = 1:6, labels = c("Dry 2018", "Wet 2018", "Dry 2019", "Wet 2019", "Dry 2020", "Wet 2020")) + labs(x = "Season", y = "Richness") + theme_bw() rich.graph # 4.2.2. SPECIES TURNOVER --------------------------------------- # Total Turnover snakec_turnover <- turnover(df = snakec, time.var = "season", species.var = "species", abundance.var = "abundance") kable(head(snakec_turnover)) # Appearance snakec_appearance <- turnover(df = snakec, time.var = "season", species.var = "species", abundance.var = "abundance", metric = "appearance") snakec_appearance # Disappearance snakec_disappearance <- turnover(df= snakec, time.var = "season", species.var = "species", abundance.var = "abundance", metric = "disappearance") snakec_disappearance # Format a compiled dataframe snakec_turnover$metric <- "total" names(snakec_turnover)[1] = "turnover" snakec_appearance$metric <- "appearance" names(snakec_appearance)[1] = "turnover" snakec_disappearance$metric <- "disappearance" names(snakec_disappearance)[1] = "turnover" snakec_allturnover <- rbind(snakec_turnover, snakec_appearance, snakec_disappearance) # Create the graph turn.graph <- ggplot(snakec_allturnover, aes(x = season, y = turnover, color = metric)) + geom_line(size = 1) + scale_x_continuous(breaks = 1:6, labels = c("Dry 2018", "Wet 2018", "Dry 2019", "Wet 2019", "Dry 2020", "Wet 2020")) + labs(x = "Season", y = "Turnover") + theme_bw() + theme(legend.position = "bottom") turn.graph # 4.2.3. MEAN RANK SHIFTS --------------------------------------- # Run the rank shift code snakec_rankshift_agg <- rank_shift(df = snakec, time.var = "season", species.var = "species", abundance.var = "abundance") snakec_rankshift_agg # Select the final time point from the returned time.var_pair snakec_rankshift_agg$season <- as.numeric(substr(snakec_rankshift_agg$year_pair, 3, 3)) # Create the graph rankshift.graph <- ggplot(snakec_rankshift_agg, aes(season, MRS)) + geom_line(size= 1) + scale_x_continuous(breaks = 1:6, labels = c("Dry 2018", "Wet 2018", "Dry 2019", "Wet 2019", "Dry 2020", "Wet 2020")) + labs(x = "Season", y = "MRS") + theme_bw() rankshift.graph # 4.2.4. RATE CHANGE -------------------------------------------- # Run the rate change code rate.res <- rate_change(snakec, time.var = "season", species.var = "species", abundance.var = "abundance") rate.res comm.res <- rate_change_interval(snakec, time.var = "season", species.var = "species", abundance.var = "abundance") comm.res quartz(8,8) rate.graph <- ggplot(comm.res, aes(interval, distance)) + geom_point() + labs(x = "Time Interval", y = "Euclidean Distance") + stat_smooth(method = "lm", se = F, size = 2) + theme_bw() rate.graph # 4.2.5. COMBINE THE PANEL GRAPHS ------------------------------- figX <- grid.arrange(rich.graph + labs(x = "Season", y = "Richness") + theme(strip.text.x = element_text(size = 14), strip.background = element_blank()) + theme(plot.margin = unit(c(0, 1, 0, 0), "cm")), turn.graph + labs(x = "Season", y = "Turnover") + theme(legend.position = "none") + theme(strip.background = element_blank(), strip.text.x = element_blank()) + scale_color_manual(values = c("#009E73", "#e79f00", "#9ad0f3")) + theme(plot.margin=unit(c(0, 1, 0, 0), "cm")), rankshift.graph + labs(x = "Season", y = "Mean Rank Shift") + theme(strip.background = element_blank(), strip.text.x = element_blank()) + theme(plot.margin = unit(c(0, 1, 0, 0), "cm")), rate.graph + labs(x = "Time interval", y = "Euclidean Distance") + theme(strip.background = element_blank(), strip.text.x = element_blank()) + theme(plot.margin = unit(c(0, 1, 0, 0), "cm")), ncol = 2) # 4.3. COMMUNITY STABILITY METRICS ----------------------------------- # 4.3.1. COMMUNITY STABILITY ------------------------------------ # Calculate community stability # média temporal dividida pelo desvio padrão temporal) snakec_stab <- community_stability(df = snakec, time.var = "season", abundance.var = "abundance") snakec_stab # 4.3.2. VARIANCE RATIO ----------------------------------------- # Calculate variance ratio, merge with stab snakec_vr <- merge(variance_ratio(df = snakec, time.var = "season", species.var = "species", abundance.var = "abundance", bootnumber = 100), snakec_stab) snakec_vr # 4.3.3. SYNCHRONY (LOREAU) ------------------------------------- # Calculate synchrony via Loreau, merge with stab # Loreau: compares the variance of aggregated species abundance with the summed variances of individual species snakec_sync_Loreau <- merge(synchrony(df = snakec, time.var = "season", species.var = "species", abundance.var = "abundance"), snakec_stab) snakec_sync_Loreau # 4.3.4. SYNCHRONY (GROSS) -------------------------------------- # Calculate synchrony via Gross, merge with stab # Gross: compares the average correlation of each individual species with the rest of the aggregated community snakec_sync_Gross <- merge(synchrony(df = snakec, time.var = "season", species.var = "species", abundance.var = "abundance", metric = "Gross"), snakec_stab) snakec_sync_Gross # 4.3.5. COMBINE THE PANEL GRAPHS ------------------------------- # Make the graphs vr_graph <-ggplot(snakec_vr, aes(x = VR, y = y)) + geom_point(size = 3) + labs(x = "Variance Ratio", y = "Community Stability") + theme_bw() + theme(text = element_text(size = 14)) vr_graph loreau_graph <- ggplot(snakec_sync_Loreau, aes(x = x, y = y)) + geom_point(size = 3) + labs(x = "Synchrony (Loreau)", y = "Community Stability") + theme_bw() + theme(text = element_text(size = 14)) loreau_graph gross_graph <- ggplot(snakec_sync_Gross, aes(x = x, y = y)) + geom_point(size = 3) + labs(x = "Synchrony (Gross)", y = "Community Stability") + theme_bw() + theme(text = element_text(size = 14)) gross_graph # Make the panel graph figY <- grid.arrange(vr_graph + labs(x = "Variance Ratio", y = "Community Stability"), loreau_graph + labs(x = "Synchrony (Loreau)", y = ""), gross_graph + labs(x = "Synchrony (Gross)", y = ""), ncol = 3)