# Platts et al 2019 Habitat availability explains variation in climate-driven range shifts across multiple # taxonomic groups. Scientific Reports. https://www.nature.com/articles/s41598-019-51582-2 library(raster) setwd("D:/REFUGE/Papers/Scientific Reports/Range margins") source("../OSGBtoXY.R") t1.starty <- 1976 # period 1: start year t1.endy <- 1990 # period 1: end year t2.starty <- 2001 # period 2: start year t2.endy <- 2015 # period 2: end year # region # ====== load("../UKHectads.rdata") region <- rasterFromXYZ(cbind(UKHectad[c("easting","northing")]+5000,1),res=10000) # temperature # =========== Tstack <- stack() for(y in t1.starty:t1.endy) { for(m in c("01","02","03","04","05","06","07","08","09","10","11","12")) { T <- raster(paste0("D:/REFUGE/Climate/UKCP09_monthly/MeanTemp/MeanTemp_",y,"-",m,"_ACTUAL.txt")) Tstack <- addLayer(Tstack,T) } } Tgrd.t1 <- mean(Tstack) Tgrd.t1 <- aggregate(Tgrd.t1,2,mean) # 10km grid Tgrd.t1 <- crop(Tgrd.t1,region) Tgrd.t1 <- mask(Tgrd.t1,region) Tval.t1 <- na.omit(values(Tgrd.t1)) rm(Tstack) par(mfrow=c(1,3)) plot(Tgrd.t1 > quantile(Tval.t1,0.5)) # hectad distributions # ==================== csvfiles <- list.files(path="../Taxonomic data/hectad data (p1p2)",pattern = "\\.csv$",full.names=TRUE) hectads <- read.csv(csvfiles[1])[c("Group","CONCEPT.NEW","NAME.NEW","Hectad","TP")] for(i in csvfiles[-1]) { if(i=="Butterflies_WSSInt1_intro_removed_updated") next hectads <- rbind( hectads, read.csv(i)[c("Group","CONCEPT.NEW","NAME.NEW","Hectad","TP")] ) } names(hectads) <- c("group","concept","latin","hectad","t") hectads$group <- as.character(hectads$group) hectads$concept <- as.character(hectads$concept) hectads$latin <- as.character(hectads$latin) hectads$hectad <- as.character(hectads$hectad) # correct name errors hectads$concept[hectads$concept=="Lep_7138"] <- "Lep_7138.ssp" hectads$concept[hectads$concept=="Lep_6508.ss"] <- "Lep_6508.ssp" hectads$concept[hectads$concept=="Lep_6457.ss"] <- "Lep_6457.ssp" hectads$concept[hectads$concept=="Lep_7972.ss"] <- "Lep_7972.ssp" hectads$concept[hectads$concept=="Hem_1656.sssp"] <- "Hem_1656.ssp" hectads$concept[hectads$concept=="Col_1223"] <- "Col_1223.f" hectads$concept[hectads$concept=="Col_7309.f"] <- "Col_1223.f" hectads$concept[hectads$concept=="Col_1224"] <- "Col_1224.f" hectads$concept[hectads$concept=="Col_7313.f"] <- "Col_1224.f" hectads$concept <- factor(hectads$concept) # remove duplicates hectads <- hectads[!duplicated(hectads),] row.names(hectads) <- seq(nrow(hectads)) # add coordinates xy <- OSGBtoXY(hectads$hectad) hectads$x <- xy$x + 5000 hectads$y <- xy$y + 5000 rm(xy) # natives # ======= hectads$native <- TRUE # identify non-natives from GBNNSIP checklist GBNNSIP <- read.csv("../Taxonomic data/GBNNSIP_with_BRC_Codes.csv") hectads$native[hectads$latin %in% GBNNSIP$NAME | hectads$concept %in% GBNNSIP$BRC_CONCEPT] <- FALSE length(unique(hectads$concept[!hectads$native])) # (finds 150 species) # identify additional 'non-natives' (and migrants, vagrants, extinct) from BRC checklist brc_native_migrant <- read.csv("../Taxonomic data/Native_Migrant_BRCStatus.csv") brc_native_migrant <- brc_native_migrant$CONCEPT.NEW[!(brc_native_migrant$NATIVE %in% c("N","NE","NNA","null"))] hectads$native[hectads$concept %in% brc_native_migrant] <- FALSE # N Native (in Mesolithic Britain or naturally colonized from Europe) # NE Native endemic # NNA Probably native # null Species wasn't in the list # exclude all others length(unique(hectads$concept[!hectads$native])) # (finds 252 additional species) # recorder effect # =============== csvfiles <- list.files(path="../Taxonomic data/2. Recorder effort",pattern = "\\.csv$",full.names=TRUE) recdata <- read.csv(csvfiles[1])[c("Group","Focal","Recorded","WellRecorded","HeavilyRecorded")] for(i in csvfiles[-1]) { recdata <- rbind( recdata, read.csv(i)[c("Group","Focal","Recorded","WellRecorded","HeavilyRecorded")] ) } recdata$Group <- as.character(recdata$Group) recdata$Focal <- as.character(recdata$Focal) # add coordinates xy <- OSGBtoXY(recdata$Focal) recdata$x <- xy$x + 5000 recdata$y <- xy$y + 5000 rm(xy) sum(is.na(recdata$Focal)) # no NAs # calculate range margins # ======================= results <- hectads[!duplicated(hectads$concept),c("group","concept","latin","native")] row.names(results) <- seq(nrow(results)) results$N.t1 <- NA # all known hectads in t1 results$N.t2 <- NA # all known hectads in t2 results$north.t1 <- NA # most northerly known hectad in t1 results$north.t2 <- NA # most northerly known hectad in t2 results$wellN.t1 <- NA # well recorded hectads in t1 results$wellN.t2 <- NA # well recorded hectads in t2 results$heavN.t1 <- NA # heavily recorded hectads in t1 results$heavN.t2 <- NA # heavily recorded hectads in t2 results$marg.t1.well <- NA # location of northern range margin for well recorded hectads in t1 results$marg.t2.well <- NA # location of northern range margin for well recorded hectads in t2 results$marg.t1N100.well <- FALSE # TRUE if at least 10 well recorded hectads to the north and to the south of t1 margin results$propwarm50.t1 <- NA # proportions of all known t1 hectads in warmest 50% of study region results$shiftkm.well <- NA results$shiftkmy.well.median <- NA results$shiftkmy.well.cilow <- NA results$shiftkmy.well.ciupp <- NA results$eligible <- FALSE marginsqs <- data.frame( matrix(nr=0,ncol=7) ) names(marginsqs) <- c("group","concept","latin","hectad","x","y","t") for(k in 1:nrow(results)) { cat(round(k/nrow(results)*100,2),"%\n") # all known hectads hectads.t1 <- hectads[hectads$concept==results$concept[k] & hectads$t==1,c("hectad","x","y")] hectads.t2 <- hectads[hectads$concept==results$concept[k] & hectads$t==2,c("hectad","x","y")] results$N.t1[k] <- nrow(hectads.t1) results$N.t2[k] <- nrow(hectads.t2) if(nrow(hectads.t1)==0 | nrow(hectads.t2)==0) next results$north.t1[k] <- max(hectads.t1$y) results$north.t2[k] <- max(hectads.t2$y) # proportion of all known hectads in warmest 50% of study region T1 <- na.omit(extract(Tgrd.t1,hectads.t1[c("x","y")])) results$propwarm50.t1[k] <- sum(T1 > quantile(Tval.t1,0.50)) / length(T1) # well recorded hectads hectads.t1.well <- hectads.t1[hectads.t1$hectad %in% recdata$Focal[recdata$Group==results$group[k] & recdata$WellRecorded==1],] hectads.t2.well <- hectads.t2[hectads.t2$hectad %in% recdata$Focal[recdata$Group==results$group[k] & recdata$WellRecorded==1],] results$wellN.t1[k] <- nrow(hectads.t1.well) results$wellN.t2[k] <- nrow(hectads.t2.well) # heavily recorded hectads hectads.t1.heav <- hectads.t1[hectads.t1$hectad %in% recdata$Focal[recdata$Group==results$group[k] & recdata$HeavilyRecorded==1],] hectads.t2.heav <- hectads.t2[hectads.t2$hectad %in% recdata$Focal[recdata$Group==results$group[k] & recdata$HeavilyRecorded==1],] results$heavN.t1[k] <- nrow(hectads.t1.heav) results$heavN.t2[k] <- nrow(hectads.t2.heav) # range margins # ~~~~~~~~~~~~~ if(results$wellN.t1[k]>=20 & results$wellN.t2[k]>=20) { # mean latitude of 10 most northern well recorded hectads top10.t1 <- sort(hectads.t1.well$y,decreasing=TRUE)[1:10] top10.t2 <- sort(hectads.t2.well$y,decreasing=TRUE)[1:10] results$marg.t1.well[k] <- mean( top10.t1 ) results$marg.t2.well[k] <- mean( top10.t2 ) # TRUE if at least 10 well recorded hectads to the north and to the south of p1 margin results$marg.t1N100.well[k] <- sum(recdata$Group==results$group[k] & recdata$WellRecorded==1 & recdata$y>=(results$marg.t1.well[k]-100000) & recdata$y=10 & sum(recdata$Group==results$group[k] & recdata$WellRecorded==1 & recdata$y<=(results$marg.t1.well[k]+100000) & recdata$y>results$marg.t1.well[k])>=10 # hectads level or north of ten most northern t1/t2 margin t1.marginsqs <- cbind( results[k,c("group","concept","latin")], hectads.t1.well[hectads.t1.well$y >= top10.t1[10],], t=1, row.names = NULL ) t2.marginsqs <- cbind( results[k,c("group","concept","latin")], hectads.t2.well[hectads.t2.well$y >= top10.t2[10],], t=2, row.names = NULL ) marginsqs <- rbind(marginsqs,t1.marginsqs,t2.marginsqs) } } # eligibility # =========== # include in study (notwithstanding habitat restrictions) results$eligible[ with(results, native & marg.t1N100.well & propwarm50.t1>=0.9 & results$north.t1<875000) ] <- TRUE # large tortoiseshell (reintroduced) results$eligible[results$latin=="Nymphalis polychloros"] <- FALSE # immigrant macromoth (missed by lookups) results$eligible[results$latin=="Trichoplusia ni"] <- FALSE # clearwings (spatial accuracy of recording affected by pheromone lures) results$eligible[results$latin %in% c("Bembecia ichneumoniformis","Sesia apiformis","Synanthedon myopaeformis","Synanthedon vespiformis")] <- FALSE # calculate range shifts # ====================== # mean shift and confidence interval around each eligible species results$shiftkm.well[results$eligible] <- (results$marg.t2.well-results$marg.t1.well)/1000 # km difference between margins nreps <- 10000 t1.samp <- 1976 + 15*runif(nreps) t2.samp <- 2001 + 15*runif(nreps) interval <- t2.samp - t1.samp for(i in 1:nrow(results)) { if(results$eligible[i]) { shift.samp <- results$shiftkm.well[i]/interval results$shiftkmy.well.median[i] <- median(shift.samp) results$shiftkmy.well.cilow[i] <- quantile(shift.samp,0.025) results$shiftkmy.well.ciupp[i] <- quantile(shift.samp,0.975) } } write.csv(results,"margins_allspp.csv",row.names=FALSE,quote=FALSE) write.csv(results[results$eligible,],"margins_eligible.csv",row.names=FALSE,quote=FALSE) write.csv(marginsqs[ marginsqs$concept %in% results$concept[results$eligible], ],"marginhectads_eligible.csv",row.names=FALSE,quote=FALSE) write.csv(hectads[ hectads$concept %in% results$concept[results$eligible], names(marginsqs)],"allhectads_eligible.csv",row.names=FALSE,quote=FALSE)