################################################################################ # August 2019 # Corresponding author: Giulia Francesca Azzurra Donati, giulia.donati@usys.ethz.ch # # Script used for figures in main text in: # Donati et al. 2019: "A process-based-model supports an association between # dispersal and the prevalence of species traits in tropical reef fish assemblages" # accepted in Ecography, August 2019 # # doi: 10.1111/ecog.04537 ################################################################################ ### Functions loading source("https://www.envidat.ch/dataset/6b06ad48-bd85-47a8-b478-0850b6def1a4/resource/eed003d7-c661-42ab-bafb-c2fd2290d94e/download/donati_et_al_dispersal_global_prevalence_traits_functions.r") #source("Donati_et_al_Dispersal_Global_Prevalence_Traits_Functions.R") # Libraries --------------------------------------------------------------- lib_vect <- c("raster","plotrix","rgeos","rgdal","sp","maptools","shape","classInt","grDevices") sapply(lib_vect,library,character.only=TRUE) # Data upload ------------------------------------------------------------- grid_pl_df<-read.csv("https://www.envidat.ch/dataset/6b06ad48-bd85-47a8-b478-0850b6def1a4/resource/2115cb5b-7557-4897-8892-f109a15eb0c1/download/grid_pl_sim_obs_trait_prop_cong.csv",h=T) #grid_pl_df<-read.csv("grid_pl_sim_obs_trait_prop_cong.csv",h=T) grid_pl_df <- grid_pl_df[,-1] spm <- read.table(url("https://www.envidat.ch/dataset/6b06ad48-bd85-47a8-b478-0850b6def1a4/resource/471350ea-cbf8-409e-a8cc-55a84fe4a8cd/download/matrix_sim_sprichness.txt/matrix_Sim_SpRichness.txt")) # lat long and species richness per simulation spm <- read.table("./matrix_Sim_SpRichness.txt") prj_coast <- CRS("+proj=cea +lat_0=0 +lon_0=0 +lat_ts=30 +a=6371228.0 +lat_1=30 +pm=-164 +units=m") newproj <- CRS("+proj=robin +lon_0=150 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +no_defs") # Modification de la grille grid<- readOGR("./", "Grid_4326") pL <- SpatialPoints(spm[,1:2],proj4string=CRS(proj4string(grid))) IDs <- over(pL, grid) IDs$Lid <- 1:dim(IDs)[1] pL <- SpatialPointsDataFrame(spm[,1:2], data=as.data.frame(IDs)) pL <- pL[which(!is.na(pL$FISHNET_ID)),] #filter only matching points ### Extract fishnet_ID FnID <- unique(pL$FISHNET_ID) FnID <- FnID[-which(FnID == 198)] grid_pl <- grid[which(grid$FISHNET_ID %in% FnID),] grid_pl <- grid_pl[-c(116,128)] coordinates <- cbind(coordinates(grid),grid$FISHNET_ID) colnames(coordinates)<-c("long","lat","FnID") ## The coastline shape file is freely available at https://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html coast <- sf::st_read("./GSHHS_l_L1_disolve.shp") ###============ Spatial congruence between simulation and observation =======### ### Data compilation for Table 2 and Fig. 1 and 2. #Spatial congruence test of overlap between simulated proportions of low and #high dispersal and observed low- and high-dispersal- related traits in assemblages (30% and 50% threshold). load("https://www.envidat.ch/dataset/6b06ad48-bd85-47a8-b478-0850b6def1a4/resource/95d1d60f-7dc9-4c51-992e-baaca842c296/download/cong_sim_obs_trait_prop.rdata") #load("cong_sim_obs_trait_prop.Rdata") ### Percentage of cells in common 30% between simulated and observed - 30% res_hot_spot <- Cons_Serv_Permut(mat_data_cong[,c(4:24)],999,seuil=30/100) Hotspot <- cbind(res_hot_spot$Mat.Hotspot, X = mat_data_cong$long,Y=mat_data_cong$lat) low_disp_sim_obs <- res_hot_spot$Resultats[which(res_hot_spot$Resultats$Var1=="psmallsim"),] high_disp_sim_obs <- res_hot_spot$Resultats[which(res_hot_spot$Resultats$Var1=="plargesim"),] res_hot_spot30 <- rbind(low_disp_sim_obs,high_disp_sim_obs) ## Percentage of cells in common 30% between simulated and observed -50% res_hot_spot2 <- Cons_Serv_Permut(mat_data_cong[,c(4:24)],999,seuil=50/100) Hotspot2 <- cbind(res_hot_spot2$Mat.Hotspot, X = mat_data_cong$long,Y=mat_data_cong$lat) low_disp_sim_obs <- res_hot_spot2$Resultats[which(res_hot_spot2$Resultats$Var1=="psmallsim"),] high_disp_sim_obs <- res_hot_spot2$Resultats[which(res_hot_spot2$Resultats$Var1=="plargesim"),] res_hot_spot50 <- rbind(low_disp_sim_obs,high_disp_sim_obs) ### Figure 1: #Global distribution maps for tropical reef fishes showing proportions of species richness for (a) low-dispersal parameter simulations (b) observed dispersal-related traits and (c) congruence zones between (a) and (b). The observed dispersal-related traits mapped are: smallest adult body sizes (left), sedentary or territorial behaviour (centre) and PLD shorter than three weeks (right). Congruence maps show overlapping zones (matching cells) between the highest proportions of low-dispersal parameter simulations and low-dispersal-related trait observations, for 30% (orange) and 50% (yellow) thresholds of hotspot congruence. ----------------------------------------------------------------- tiff(filename = "./Figure_2",width=21,height=18,units="cm",res=300,compression="lzw",type="cairo") layout(mat=rbind(c(1,1,1,1,1,1,1,1,2),c(1,1,1,1,1,1,1,1,2),c(1,1,1,1,1,1,1,1,2), c(1,1,1,1,1,1,1,1,2),c(3,3,3,4,4,4,5,5,5),c(3,3,3,4,4,4,5,5,5), c(6,6,6,7,7,7,8,8,8),c(9,9,9,10,10,10,11,11,11),c(9,9,9,10,10,10,11,11,11), c(12,12,12,13,13,13,14,14,14))) col_vec <- colorRampPalette(c("deepskyblue", "darkorchid2","firebrick1"))(6) par(mar=c(2,2,2.2,1))# (a) Plot_the_world(var=na.omit(grid_pl_df$psmallsim), coast=coast,grid=grid,col_vec=col_vec, cex_axy=0.9,cex_axx=0.9,xlim=c(-165,167),ylim=c(-50,50)) mtext(text="a) Simulations",at=60,side=2,las=2,line=-5,cex=0.8) par(mar =c(2,3,4,3)) # legend (a) Plot_lengend_bouy(var=na.omit(grid_pl_df$psmallsim),color=col_vec, title="% lowest\n dispersing",h=FALSE,cex.main=1,cex.axis=.8) par(mar=c(1,2,2,0)) # (b) Plot_the_world(var=grid_pl_df$psmallobs10,coast=coast,grid=grid,col_vec=col_vec, axey=TRUE,axex=TRUE,xlim=c(-168,177),ylim= c(-55,55),lwdgrid=0.3, labey=c(-45,-15,15,45)) mtext(text = " b) Observations",at=75,side=2,las=2,line=-5.8,cex=0.8) par(mar=c(1,1,2,1)) Plot_the_world(var=grid_pl_df$psedobs,coast=coast,grid=grid,col_vec=col_vec,axey=FALSE, axex=TRUE,xlim=c(-168,177),ylim= c(-55,55),lwdgrid=0.3,labey=c(-45,-15,15,45)) par(mar=c(1,0,2,2)) Plot_the_world(var = grid_pl_df$pshort_pld,coast=coast,grid=grid,col_vec=col_vec,axey=FALSE, axex=TRUE,xlim=c(-168,177),ylim=c(-55,55),lwdgrid=0.3,labey=c(-45,-15,15,45)) par(mar =c(2,3,2.5,1))# legends (b) Plot_lengend_bouy(var=grid_pl_df$psmallobs10,color=col_vec,title="% smallest-bodied", zlevels=5,cex.axis=.6,cex.main=1) par(mar =c(2,3,2.5,2)) Plot_lengend_bouy(var=grid_pl_df$psedobs,color=col_vec,title="% sedentary or territorial", zlevels=5,cex.axis=.6,cex.main=1) par(mar =c(2,0,2.5,2)) Plot_lengend_bouy(var=grid_pl_df$pshort_pld,color=col_vec,title="% PLD (median) < 3 weeks", zlevels=5,cex.axis=.6,cex.main=1) par(mar=c(1,2,2,0))# (c) carto_congruence_2(data="Hotspot", data2="Hotspot2", x="psmallsim",index2="psmallobs10",names_fig="a", xlim=c(-168,177),ylim= c(-55,55),cex_point=0.4,cex_text_fig=1,cex_legend=0.5, cex_names_fig=0.1,ncol=2, bathy=FALSE,coast=coast,labey=c(-45,-15,15,45), col=c("darkorange1","gold")) mtext(text = " c) Congruence zones",at=75,side=2,las=2,line=-8.2,cex=0.8) par(mar=c(1,1,2,1)) carto_congruence_2(data="Hotspot", data2="Hotspot2", x="psmallsim",index2="psedobs",names_fig="a", xlim=c(-168,177),ylim= c(-55,55),cex_point=0.7,cex_text_fig=1,cex_legend=0.5, cex_names_fig=0.1, ncol=2, bathy=FALSE,coast=coast,axey=FALSE) par(mar=c(1,0,2,2)) carto_congruence_2(data="Hotspot", data2="Hotspot2", x="psmallsim",index2="pshort_pld",names_fig="a", xlim=c(-168,177),ylim= c(-55,55),cex_point=0.7,cex_text_fig=1,cex_legend=0.5,cex_names_fig=0.1, ncol=2, bathy=FALSE,coast=coast,axey=FALSE) # legends (c) par(mar=c(0.75,0,1,2)) emptyplot(ylim=c(0,1),xlim=c(0,1)) legend("center",leg=c("30%","50%"),pch=15,col=c("darkorange1","gold"),bg="white", pt.cex=0.9,cex=1,lwd=0.2,box.lwd=0.6,lty=0,ncol=1,bty="n",x.intersp=0.5, title="low dispersal and smallest-bodied") par(mar=c(0.75,0,1,2)) emptyplot(ylim=c(0,1),xlim=c(0,1)) legend("center",leg=c("30%","50%"),pch=15,col=c("darkorange1","gold"),bg="white", pt.cex=0.9,cex=1,lwd=0.2,box.lwd=0.6,lty=0,ncol=1,bty="n",x.intersp=0.5, title="low dispersal and restricted home range") par(mar=c(0.75,0,1,2)) emptyplot(ylim=c(0,1),xlim=c(0,1)) legend("center",leg=c("30%","50%"),pch=15,col=c("darkorange1","gold"),bg="white", pt.cex=0.9,cex=1,lwd=0.2,box.lwd=0.6,lty=0,ncol=1,bty="n",x.intersp=0.5, title="low dispersal and PLD (median) < 3 weeks") dev.off() # Figure 2: Global distribution maps for tropical reef shes showing proportions of species richness for (a) high-dispersal parameter simula- tions (b) observed dispersal-related traits and (c) congruence zones between (a) and (b). e observed dispersal-related traits mapped are: the largest adult-body sizes (left), highly mobile behaviour (centre), and PLD longer than six weeks (right). Congruence maps show overlap- ping zones (matching cells) between the highest proportions of high-dispersal parameter simulations and high-dispersal-related trait obser- vations, for 30% (orange) and 50% (yellow) thresholds of hotspot congruence.----------------------------------------------------------------- tiff(filename = "./Figure_2",width =21, height = 18, units = "cm", res=300,compression = "lzw",type="cairo") layout(mat=rbind(c(1,1,1,1,1,1,1,1,2),c(1,1,1,1,1,1,1,1,2),c(1,1,1,1,1,1,1,1,2), c(1,1,1,1,1,1,1,1,2),c(3,3,3,4,4,4,5,5,5),c(3,3,3,4,4,4,5,5,5), c(6,6,6,7,7,7,8,8,8),c(9,9,9,10,10,10,11,11,11),c(9,9,9,10,10,10,11,11,11), c(12,12,12,13,13,13,14,14,14))) par(mar=c(2,2,2.2,1)) col_vec <- colorRampPalette(c("deepskyblue", "darkorchid2","firebrick1"))(6) # (a) Plot_the_world(var=na.omit(grid_pl_df$plargesim), coast=coast,grid=grid,col_vec=col_vec, cex_axy=0.9,cex_axx = 0.9,xlim=c(-165,167),ylim= c(-50,50)) mtext(text = "a) Simulations",at=60,side=2,las=2,line=-5,cex=0.8) par(mar =c(2,3,4,3))# legend (a) Plot_lengend_bouy(var=na.omit(grid_pl_df$plargesim),color=col_vec,title="% highest\n dispersing", h=FALSE,cex.main=1,cex.axis=.8) par(mar=c(1,2,2,0)) # (b) Plot_the_world(var=grid_pl_df$plargeobs90,coast=coast,grid=grid,col_vec=col_vec, axey=TRUE,axex=TRUE,xlim=c(-168,177),ylim= c(-55,55),lwdgrid=0.3, labey=c(-45,-15,15,45)) mtext(text=" b) Observations",at=75,side=2,las=2,line=-5.8,cex=0.8) par(mar=c(1,1,2,1)) Plot_the_world(var=grid_pl_df$pmobobs,coast=coast,grid=grid,col_vec=col_vec, axey=FALSE,axex=TRUE,xlim=c(-168,177),ylim= c(-55,55),lwdgrid=0.3, labey=c(-45,-15,15,45)) par(mar=c(1,0,2,2)) Plot_the_world(var=grid_pl_df$plong_pld,coast=coast,grid=grid,col_vec=col_vec, axey=FALSE,axex=TRUE,xlim=c(-168,177),ylim= c(-55,55),lwdgrid=0.3, labey=c(-45,-15,15,45)) par(mar =c(2,3,2.5,1))# legends (b) Plot_lengend_bouy(var=grid_pl_df$plargeobs90,color=col_vec,title="% largest-bodied", zlevels=5,cex.axis=.6,cex.main=1) par(mar =c(2,3,2.5,2)) Plot_lengend_bouy(var=grid_pl_df$pmobobs,color=col_vec,title="% mobile-highly mobile between reefs", zlevels=5,cex.axis=.6,cex.main=1) par(mar =c(2,0,2.5,2)) Plot_lengend_bouy(var=grid_pl_df$plong_pld,color=col_vec,title="% PLD (median) > 6 weeks", zlevels=5,cex.axis=.6,cex.main=1) par(mar=c(1,2,2,0))# (c) carto_congruence_2(data="Hotspot",data2="Hotspot2",x="plargesim",index2="plargeobs90",names_fig="a", xlim=c(-168,177),ylim= c(-55,55),cex_point=0.4,cex_text_fig=1,cex_legend=0.5, cex_names_fig=0.1,ncol=2, bathy=FALSE,coast=coast,labey=c(-45,-15,15,45)) mtext(text=" c) Congruence zones",at=75,side=2,las=2,line=-8.2,cex=0.8) par(mar=c(1,1,2,1)) carto_congruence_2(data="Hotspot",data2="Hotspot2",x="plargesim",index2="pmobobs",names_fig="a", xlim=c(-168,177),ylim= c(-55,55),cex_point=0.7,cex_text_fig=1,cex_legend=0.5, cex_names_fig=0.1,ncol=2,bathy=FALSE,coast=coast,axey=FALSE) par(mar=c(1,0,2,2)) carto_congruence_2(data="Hotspot",data2="Hotspot2",x="plargesim",index2="plong_pld",names_fig="a", xlim=c(-168,177),ylim=c(-55,55),cex_point=0.7,cex_text_fig=1, cex_legend=0.5,cex_names_fig=0.1,ncol=2,bathy=FALSE,coast=coast,axey=FALSE) par(mar=c(0.75,0,1,2)) # legends (c) emptyplot(ylim=c(0,1),xlim=c(0,1)) legend("center",leg=c("30%","50%"),pch=15,col=c("darkorange1","gold"),bg="white", pt.cex=0.9,cex=1,lwd=0.2,box.lwd=0.6,lty=0,title="high dispersal and large-bodied", ncol=1,bty="n",x.intersp=0.5) par(mar=c(0.75,0,1,2)) emptyplot(ylim=c(0,1),xlim=c(0,1)) legend("center",leg=c("30%","50%"),pch=15,col=c("darkorange1","gold"),bg="white", pt.cex=0.9,cex=1,lwd=0.2,box.lwd=0.6,lty=0,title="high dispersal and wide home range", ncol=1,bty="n",x.intersp=0.5) par(mar=c(0.75,0,1,2)) emptyplot(ylim=c(0,1),xlim=c(0,1)) legend("center",leg=c("30%","50%"),pch=15,col=c("darkorange1","gold"),bg="white", pt.cex=0.9,cex=1,lwd=0.2,box.lwd=0.6,lty=0,title="high dispersal and PLD (median) > 6 weeks", ncol=1,bty="n",x.intersp=0.5) dev.off() # Figure 3: Frequency of distribution of the dispersal parameter considered in the simulations (a) and observed frequency distribution for adult body size (b), median PLD (c), and home range category (d). The number of species is biased toward more smaller species, with lower median PLD and of sedentary home range, which corresponds to the higher frequency of simulated species with lower dispersal.--------------------------------------------------------------- spt_gaspar<-read.csv("./sp_ecological_trait_data.csv")#traitdata tiff(filename="./Figure_3.tiff",width=20,height=17,units="cm",res=300,compression="lzw",type="cairo") par(mfrow=c(2,2),lwd=0.5) par(mar=c(5.1,4,2.2,1)) hist(spt$x,50,main="",xlab="Dispersal parameter (d)",col="grey",cex.main=0.8,cex.lab=0.8,cex.axis=0.8,lwd=0.5,las=1) mtext("a)",side=3,cex=.8,font=1,las=1,line=1,at=1) hist(spt_gaspar$Max_length,main="",70,xlab="Adult body size (cm)",ylab="",xlim=c(0,200), ylim=c(0,1500), cex.main=0.8,col="grey",cex.lab=0.8,cex.axis=0.8,lwd=0.5,las=1) mtext("b)",side=3,cex=.8,font=1,las=1,line=1,at=0) hist(spt_gaspar$vect_mPLD, main="",10,xlab="PLD median (days)",xlim=c(0,100),ylim=c(0,1500), cex.main=0.8,col="grey",cex.lab=0.8,cex.axis=0.8,lwd=0.5,las=1) mtext("c)",side=3,cex=.8,font=1,las=1,line=1,at=0) spt_gaspar$Home_Range <- gsub("Sed","sed",as.character(spt_gaspar$Home_Range)) spt_gaspar$Home_Range <- gsub("VMob","h-mob",as.character(spt_gaspar$Home_Range))# order is important spt_gaspar$Home_Range <- gsub("Mob","mob",as.character(spt_gaspar$Home_Range)) spt_gaspar$Home_Range <- ordered(spt_gaspar$Home_Range,levels=c("sed","mob","h-mob")) counts <- table(spt_gaspar$Home_Range) barplot(counts,main="",xlab="Home range categories",cex.main=0.8,col="grey",cex.lab=0.8, cex.axis=0.8,lwd=0.5,las=1,width=0.1,space=0.4) mtext("d)",side=3,cex=.8,font=1,las=1,line=1,at=0) dev.off()