The script extracts climatic data for populations in the experimental site (clonal bank: CB) and their natural locations. Climatic data from ClimateDT are available in two formats: point-based layers and raster data.
For the populations’ origin areas, we used climatic data from the
reference period (1901–1950) in raster format. For the CB in Spain, we
used specific time periods corresponding to each phenotypic trait
measurement:
- Shoot volume, stem length, spring elongation, total open mature
strobili (male): 1992–2012
- Leaf thickness: 1992–2021
Climatic data from ClimateDT are available as point-based layers and raster data. For the CB location, we used point-based data to obtain climate information over the period between the CB’s establishment and the year of trait measurement. Raster data, which offer more climate models for future projections, were applied to analyse the range-wide dataset. This script aims to verify consistency between point-based and raster data for CB populations, to determine whether raster data can be used to validate models built on point-based data.
For each trait, the climatic period corresponds to the years from the
establishment of the clonal bank (1992) to the latest year of trait
measurement. One key difference between raster and point-based data lies
in their formats: - Point-based data: one value per year for each
variable at each population/pixel selected.
- Raster data: less precise, representing average values over the
selected period.
Bioclimatic (bio) and other climate variables are derived from monthly minimum temperature (tmin), maximum temperature (tmax), and precipitation. To obtain an accurate mean for bio variables over a period (e.g., bio1 from 1901–1950), we first average tmin, tmax, and precipitation for each month individually, resulting in one value per month (January, February, etc.) over the period of interest.
Using these monthly averages, we apply the biovars function from the Dismo R package to compute the 19 bioclimatic variables for the specified period (e.g., 1901–1950). For other climate variables, such as AHM, over the same period, we follow ClimateDT’s protocol by incorporating the newly generated bio1, bio10, bio12, etc., values from biovars.
meta_data <- read.csv("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Populations/Populations_CG_regions.csv",h=T,sep=";",dec=",")
First, we loaded the tmin, tmax, and precipitation data from the establishment of the clonal bank to the last measurement (1992–2021), extracted from ClimateDT.
Pre_past_climatic_data_point <- read.csv("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Climatic_data/Common_garden_pop/clim_data_pop_CG.csv",h=T,sep=";",dec=",")
Calculation of the bioclimatic variables using biovars:
climatic_data_filtered <- data.frame(Pre_past_climatic_data_point) %>%
mutate(across(6:41, as.numeric))#Pass the variables in numeric
#Do the mean of each year of data for each climatic variables and that's for each ID (pop)
climatic_data_filtered_mean <- climatic_data_filtered %>%
group_by(ID) %>%
summarize(across((5):(40), ~ mean(.))) ##-1 because, ID is no longer a columns but a column names so it does not count it as a column
#We need to extract the ID, longi and lati from the dataframe with all the rows
climatic_data_mean <- climatic_data_filtered %>%
group_by(ID) %>% #Groups the row by ID
slice_head(n = 1) %>% #Keep only the first row of each groups
dplyr::select(c(1:3))#Select only the 3 first columns
#Add the longitude and latitude variables
Mean_pre_variables <- data.frame(climatic_data_mean[,c(1:3)],climatic_data_filtered_mean[,-1]) #-1 because ID is already in
data_climatic <- data.frame(Mean_pre_variables)
#We extracted the prec, tmn and tmx from the inpute data
prec <- data.frame(t(data_climatic %>% dplyr::select(starts_with("prc"))))
tmn <- data.frame(t(data_climatic %>% dplyr::select(starts_with("tmn"))))
tmx <- data.frame(t(data_climatic %>% dplyr::select(starts_with("tmx"))))
##Bio calculation
Bio_data <- data.frame(biovars(prec,tmn,tmx))
Raster data were obtained from Maurizio Marchi (IBBR-CNR), one of the creators of ClimateDT, because raster data are not directly available on ClimateDT.
#Load past and future raster data
climatic_variables_to_keep <-c("bio1", "bio12", "bio15","bio2","bio4","bio9") #Lubset of retained clim var, see script processing clim data
list_raster <- c("past_climatic_data_raster")
for(x in 1:length(list_raster)){
name_clim <- c("past")#Name final raster
var <- list_raster[x]#Name of the raster
load(paste0("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Genomic_offset/RDA/",var,".Rdata")) #Load data
raster_clim <- get(var)[[climatic_variables_to_keep]]#Keep only var in clim_variables to keep
names(raster_clim)=c("Annual_Tc","Annual_P","P_Seasonality","Diurnal_range_Tc","Tc_Seasonality","Tc_driest_quarter") #Change the name to match the names of the scaled and center values
assign(paste0("raster_",name_clim,"_clim"),raster_clim)
}
#First we need to extract the climatic data
list_coord <- c("coords_pop")
list_name <- c("natural_pop")
coords_pop <- data.frame(apply(meta_data[,c(3:4)], 2, as.numeric))#We need to have longitude then latitude
coords_CG <- data.frame(meta_data[,c(1,2)],Longitude=c("-4.0125"),Latitude=c("40.91056")) #Coordinates of the clonal bank
coords_CG_f <- data.frame(coords_CG[,-c(1,2)])
coords_CG_f$Longitude<- as.numeric(coords_CG_f$Longitude);coords_CG_f$Latitude<- as.numeric(coords_CG_f$Latitude)
for(x in 1: length(list_name)){
coord <- get(list_coord[x])
name <- list_name[x]
clim_data <- raster::extract(raster_past_clim, coord)
clim_data_order <- data.frame(clim_data[,c(1,4,5,6,2,3)])
assign(paste0("clim_df_",name),clim_data_order)
}
#If this is the period 1901_1950 then
layer_data <- data.frame(climatic_data_filtered_mean$ID,meta_data$Region,Bio_data[,c(1,2,4,9,12,15)])
colnames(layer_data) <- c("Population","Region","bio1","bio2","bio4","bio9","bio12","bio15")
colnames(layer_data)[3:8] <- paste0(colnames(layer_data)[3:8], "_layer_of_points")
#Raster data for period 1901-1950
raster_data <- data.frame(climatic_data_filtered_mean$ID,meta_data$Region,clim_df_natural_pop)
colnames(raster_data) <- c("Population","Region","bio1","bio2","bio4","bio9","bio12","bio15")
colnames(raster_data)[3:8] <- paste0(colnames(raster_data)[3:8], "_raster")
#General data frame with the 2 methods
data_both_method <- data.frame(merge(layer_data,raster_data,"Population"))
bioclim_vars <- c("bio1","bio2","bio4","bio9","bio12","bio15")
for (var in bioclim_vars) {
#Extract the data for the bioclimatic variable for the two periods
data_subset <- data.frame(data_both_method[, paste0(var, "_layer_of_points")], data_both_method[, paste0(var,"_raster")]) #For each variable in bioclim_vars, we extract the two periods and created a data.frame with only these two variables.
colnames(data_subset) <- c("var1", "var2")
#Create scatterplot to compare the bioclimatic variable between the two periods
scatterplot_plot <-
ggplot(data_subset,aes(x=data_subset[,1],y=data_subset[,2])) +
geom_point(aes(color=data_both_method$Region.x),size=3)+
geom_abline(intercept = 0, slope = 1, color="gray60")+
scale_colour_manual(name = "Regions",
values = c("orangered3","gold2","darkorchid3","navyblue","turquoise2","darkgreen")) +
xlab("Layer of points")+
ylab("Raster")+
ggtitle(paste0("Comparison climatic values ",var," methods for the
26 Taxus baccata populations planted in the clonal bank "))+
theme_set(theme_bw())
plot(scatterplot_plot)
#Save the plots
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/climate/Common_garden_comparison_layer_raster/Screeplot_comparison_",var,"_raster_layer.pdf"));print(scatterplot_plot);dev.off()
}
Overall, we observed that for the majority of populations, values from raster and point-based layers are nearly identical. These results suggest that the two data types could be used interchangeably.
meta_data_wo_regions <- read.csv("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Populations/Populations_common_garden.csv",h=T,sep=";",dec=",")
#Load past and future raster data
climatic_variables_to_keep <-c("bio1", "bio12", "bio15","bio2","bio4","bio9") #Subset of retained clim var, see script processing clim data
list_raster <- c("past_climatic_data_raster")
for(x in 1:length(list_raster)){
name_clim <- c("present")#Name final raster
var <- list_raster[x]#Name of the raster
load(paste0("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Genomic_offset/RDA/",var,".Rdata")) #Load data
raster_clim <- get(var)[[climatic_variables_to_keep]]#Keep only var in clim_variables to keep
names(raster_clim)=c("Annual_Tc","Annual_P","P_Seasonality","Diurnal_range_Tc","Tc_Seasonality","Tc_driest_quarter") #Change the name to match the names of the scaled and center values
assign(paste0("raster_",name_clim,"_clim"),raster_clim)
}
#First we need to extract the climatic data
list_coord <- c("coords_pop")
list_name <- c("natural_pop")
#To ensure consistency, we standardised by using the same values as the climatic data used in the GEA models
load("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Climatic_data/new_selection/scale_env_value_new_cli.Rdata")
load("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Climatic_data/new_selection/center_env_value_new_cli.Rdata")
coords_pop <- data.frame(apply(meta_data_wo_regions[,c(2:3)], 2, as.numeric))#We need to have longitude then latitude
for(x in 1: length(list_name)){
coord <- get(list_coord[x])
name <- list_name[x]
clim_data <- raster::extract(raster_present_clim, coord)
clim_data_order <- data.frame(clim_data[,c(1,4,5,6,2,3)])#Put it in the same order as center_env and scale object
#Standardisation of the environmental variables
clim_data_scale <- as.data.frame(scale(clim_data_order, center=center_env_value_new_cli, scale=scale_env_value_new_cli))
assign(paste0("clim_df_",name),clim_data_scale)
}
save(clim_df_natural_pop,file="C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/validation_GO/climatic_data/clim_df_natural_pop.Rdata")
We can load the data and separate it by periods:
Climatic_data_CG <- read.csv("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Climatic_data/Common_garden_pop/CG_climatic_data_T.csv",h=T,sep=";",dec=",")
clim_CG_2012 <- Climatic_data_CG[Climatic_data_CG$Year <= 2012,]
clim_CG_2021 <- Climatic_data_CG
We can calculate the bioclim variables with the biovar function:
list_clim_df <- c("clim_CG_2012","clim_CG_2021")
for(x in 1: length(list_clim_df)){
clim_df <- get(list_clim_df[x])
name <- list_clim_df[x]
climatic_data_filtered <- data.frame(clim_df) %>%
mutate(across(6:41, as.numeric))#Pass the variables in numeric
#do the mean of each year of data for each climatic variables and that's for each ID (pop)
climatic_data_filtered_mean <- climatic_data_filtered %>%
group_by(ID) %>%
summarize(across((5):(40), ~ mean(.))) ##-1 because, ID is no longer a columns but a column names so it does not count it as a column
#we need to extract the ID, longi and lati from the dataframe with all the rows
climatic_data_mean <- climatic_data_filtered %>%
group_by(ID) %>% #Groups the row by ID
slice_head(n = 1) %>% #Keep only the first row of each groups
dplyr::select(c(1:3))#Select only the 3 first columns
#Add the longitude and latitude variables
Mean_pre_variables <- data.frame(climatic_data_mean[,c(1:3)],climatic_data_filtered_mean[,-1]) #-1 because ID is already in
data_climatic <- data.frame(Mean_pre_variables)
#We extracted the prec, tmn and tmx from the inpute data
prec <- data.frame(t(data_climatic %>% dplyr::select(starts_with("prc"))))
tmn <- data.frame(t(data_climatic %>% dplyr::select(starts_with("tmn"))))
tmx <- data.frame(t(data_climatic %>% dplyr::select(starts_with("tmx"))))
##Bio calculation
Bio_data <- data.frame(biovars(prec,tmn,tmx))
assign(paste0("df_",name),Bio_data)
}
Scaling the values of the retained bio variables:
list_clim_name <- c("CG_2012","CG_2021")
for(x in 1: length(list_clim_df)){
name <- list_clim_name[x]
df <- get(paste0("df_clim_",name))
df_retained <- data.frame(df[,c(1,2,4,9,12,15)])
colnames(df_retained) <- row.names(data.frame(scale_env_value_new_cli))
assign(paste0("clim_df_",name),df_retained)
#scale
df_scale <- as.data.frame(scale(df_retained, center=center_env_value_new_cli, scale=scale_env_value_new_cli))
assign(paste0("clim_df_scale_",name),df_scale)
}
We can analyse the climate of the clonal bank (CB) in comparison to that of the populations to determine whether the CB represents an extreme environment relative to the populations. This investigation will help us understand climatic differences and assess how they may influence the growth and phenotypic traits of populations in the common garden.
#CB climatic data 1992-
load("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/validation_GO/climatic_data/clim_df_CG_2012.Rdata")
load("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/validation_GO/climatic_data/clim_df_CG_2021.Rdata")
clim_df_CG_2012_tot <- data.frame(Population=c("CG_2012"),Region="CG_2012",clim_df_CG_2012)
clim_df_CG_2021_tot <- data.frame(Population=c("CG_2021"),Region="CG_2021",clim_df_CG_2021)
#Pop climatic data 1901-1950
colnames(raster_data) <- colnames(clim_df_CG_2012_tot)
data_clim_tot <- rbind(clim_df_CG_2012_tot,clim_df_CG_2021_tot,raster_data)
#Do the violin plot in a loop for each climatic variables
#Input data for the loop
data <- data_clim_tot
regions <-data_clim_tot$Region
bioclim_vars <- colnames(clim_df_CG_2012_tot[,-c(1,2)])
unit_var <- c("°c","°c","°c","°c","mm","mm")
#Results: the loop
for (x in 1:length(bioclim_vars)) {
var <- bioclim_vars[x]
unit <- unit_var[x]
data_subset <- data.frame(data[, paste0(var)]) #For each variable in bioclim_vars, we extract the two periods and created a data.frame with only these two variables.
colnames(data_subset) <- c("var1")
#Create scatterplot to compare the bioclimatic variable between the two periods
plot <- ggplot(data_subset, aes(x =var,y= var1)) +#Var is the values and var1 is just the name of the variable
geom_violin(trim = FALSE) +
scale_colour_manual(name = "Region",
values = c("orangered3","gold2","darkorchid3","blue","turquoise2","red","black","darkgreen"))+
geom_point(aes(color = regions), shape = 16,size=3, position = position_jitter(seed = 1, width = 0.2)) +
labs(colour="Regions",x=var,y=paste0("values ","(",unit,")"),title=paste0("Violin Plot of ",var," differences across populations"))+
#theme(plot.title = element_text(hjust = 0.5))
theme_classic() + #Removes the grey background
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12),
legend.key = element_rect(fill = "transparent", color = NA),
panel.grid.major = element_line(color = "grey"), # Grille principale noire
panel.grid.minor = element_line(color = "grey"), # Grille secondaire noire
panel.grid.major.x = element_line(color = "grey"), # Grille verticale noire
panel.grid.major.y = element_line(color = "grey") # Grille horizontale noire
)
#Add any additional customization or saving the plot if needed
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/climate/Common_garden _climate/violin_plot/violin_plot_",var,".pdf"));print(plot)
;dev.off()
#Print the plot
print(plot)
}
Interpretation: The climatic range among populations is substantial, particularly in terms of annual precipitation and temperature. For the clonal bank (CB), we observed that, compared to other populations, the annual temperature is relatively high, especially during the driest quarter. Temperature seasonality and diurnal temperature range are also notable.
Regarding precipitation, annual values are relatively low, while seasonal precipitation generally follows the mean of the other populations.
These results suggest that the climate of the CB could impose higher aridity constraints and greater temperature seasonality due to high annual temperature variability. This implies that adaptation to drought resistance and later growth or reproductive phenology may be advantageous in the CB climate, given the combination of high temperatures during the driest period, low precipitation, and significant seasonal temperature fluctuations, which may also result in late frost events.
devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.3.2 (2023-10-31 ucrt)
## os Windows 11 x64 (build 26100)
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate French_France.utf8
## ctype French_France.utf8
## tz Europe/Paris
## date 2025-09-22
## pandoc 3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## bslib 0.6.1 2023-11-28 [1] CRAN (R 4.3.2)
## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.3.1)
## cli 3.6.1 2023-03-23 [1] CRAN (R 4.3.1)
## codetools 0.2-19 2023-02-01 [2] CRAN (R 4.3.2)
## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.1)
## corrplot * 0.92 2021-11-18 [1] CRAN (R 4.3.2)
## devtools 2.4.5 2022-10-11 [1] CRAN (R 4.3.3)
## digest 0.6.33 2023-07-07 [1] CRAN (R 4.3.1)
## dismo * 1.3-14 2023-05-21 [1] CRAN (R 4.3.2)
## dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.3.2)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.1)
## evaluate 0.23 2023-11-01 [1] CRAN (R 4.3.2)
## farver 2.1.2 2024-05-13 [1] CRAN (R 4.3.3)
## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.1)
## forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.3.2)
## fs 1.6.3 2023-07-20 [1] CRAN (R 4.3.1)
## generics 0.1.4 2025-05-09 [1] CRAN (R 4.3.2)
## ggplot2 * 3.5.2 2025-04-09 [1] CRAN (R 4.3.2)
## glue 1.7.0 2024-01-09 [1] CRAN (R 4.3.2)
## gtable 0.3.6 2024-10-25 [1] CRAN (R 4.3.3)
## highr 0.10 2022-12-22 [1] CRAN (R 4.3.1)
## hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.2)
## htmltools 0.5.7 2023-11-03 [1] CRAN (R 4.3.2)
## htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.3.2)
## httpuv 1.6.13 2023-12-06 [1] CRAN (R 4.3.2)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.3.1)
## jsonlite 1.8.8 2023-12-04 [1] CRAN (R 4.3.2)
## knitr 1.45 2023-10-30 [1] CRAN (R 4.3.2)
## labeling 0.4.3 2023-08-29 [1] CRAN (R 4.3.1)
## later 1.3.2 2023-12-06 [1] CRAN (R 4.3.2)
## lattice 0.22-5 2023-10-24 [2] CRAN (R 4.3.2)
## lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.2)
## lubridate * 1.9.3 2023-09-27 [1] CRAN (R 4.3.2)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.1)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.1)
## mime 0.12 2021-09-28 [1] CRAN (R 4.3.1)
## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.3.2)
## munsell 0.5.1 2024-04-01 [1] CRAN (R 4.3.3)
## pillar 1.10.2 2025-04-05 [1] CRAN (R 4.3.3)
## pkgbuild 1.4.3 2023-12-10 [1] CRAN (R 4.3.2)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.2)
## pkgload 1.3.4 2024-01-16 [1] CRAN (R 4.3.2)
## profvis 0.3.8 2023-05-02 [1] CRAN (R 4.3.2)
## promises 1.2.1 2023-08-10 [1] CRAN (R 4.3.2)
## purrr * 1.0.2 2023-08-10 [1] CRAN (R 4.3.2)
## R6 2.6.1 2025-02-15 [1] CRAN (R 4.3.3)
## raster * 3.6-26 2023-10-14 [1] CRAN (R 4.3.2)
## Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.3.1)
## readr * 2.1.5 2024-01-10 [1] CRAN (R 4.3.2)
## remotes 2.5.0 2024-03-17 [1] CRAN (R 4.3.3)
## rlang 1.1.3 2024-01-10 [1] CRAN (R 4.3.2)
## rmarkdown 2.25 2023-09-18 [1] CRAN (R 4.3.1)
## rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.2)
## sass 0.4.8 2023-12-06 [1] CRAN (R 4.3.2)
## scales 1.3.0 2023-11-28 [1] CRAN (R 4.3.2)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.2)
## shiny 1.8.0 2023-11-17 [1] CRAN (R 4.3.2)
## sp * 2.1-3 2024-01-30 [1] CRAN (R 4.3.2)
## stringi 1.8.3 2023-12-11 [1] CRAN (R 4.3.2)
## stringr * 1.5.1 2023-11-14 [1] CRAN (R 4.3.2)
## terra 1.7-71 2024-01-31 [1] CRAN (R 4.3.2)
## tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.3.2)
## tidyr * 1.3.1 2024-01-24 [1] CRAN (R 4.3.2)
## tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.3.3)
## tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.3.2)
## timechange 0.3.0 2024-01-18 [1] CRAN (R 4.3.2)
## tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.3.2)
## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.3.2)
## usethis 2.2.3 2024-02-19 [1] CRAN (R 4.3.2)
## vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.2)
## withr 3.0.2 2024-10-28 [1] CRAN (R 4.3.3)
## xfun 0.41 2023-11-01 [1] CRAN (R 4.3.2)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.2)
## yaml 2.3.8 2023-12-11 [1] CRAN (R 4.3.2)
##
## [1] C:/Users/tfrancisco/AppData/Local/R/win-library/4.3
## [2] C:/Users/tfrancisco/AppData/Local/Programs/R/R-4.3.2/library
##
## ──────────────────────────────────────────────────────────────────────────────