sessionInfo()
## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8 
## [2] LC_CTYPE=English_United Kingdom.utf8   
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.utf8    
## 
## time zone: Europe/Berlin
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] rnaturalearthdata_1.0.0 rnaturalearth_1.2.0     ggforce_0.5.0          
##  [4] sf_1.0-24               usdm_2.1-7              geosphere_1.6-8        
##  [7] geodata_0.6-9           terra_1.8-93            dismo_1.3-16           
## [10] raster_3.6-32           sp_2.2-1                ecodist_2.1.3          
## [13] ggrepel_0.9.6           readr_2.1.6             vegan_2.7-2            
## [16] permute_0.9-10          poppr_2.9.8             dartR_2.9.9.5          
## [19] dartR.data_1.2.2        dplyr_1.2.1             ggplot2_4.0.2          
## [22] adegenet_2.1.11         ade4_1.7-23             vcfR_1.15.0            
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3  rstudioapi_0.18.0   jsonlite_2.0.0     
##   [4] magrittr_2.0.3      farver_2.1.2        rmarkdown_2.30     
##   [7] fields_17.3         vctrs_0.7.3         SNPRelate_1.42.0   
##  [10] mmod_1.3.3          htmltools_0.5.8.1   sass_0.4.10        
##  [13] KernSmooth_2.23-26  bslib_0.10.0        gsubfn_0.7         
##  [16] plyr_1.8.9          cachem_1.1.0        igraph_2.2.2       
##  [19] mime_0.13           lifecycle_1.0.5     iterators_1.0.14   
##  [22] pkgconfig_2.0.3     gap_1.15.2          Matrix_1.7-4       
##  [25] R6_2.6.1            fastmap_1.2.0       rbibutils_2.4.1    
##  [28] shiny_1.12.1        digest_0.6.37       polysat_1.7-7      
##  [31] GGally_2.4.0        patchwork_1.3.2     gdistance_1.6.5    
##  [34] gdata_3.0.1         polyclip_1.10-7     mgcv_1.9-4         
##  [37] compiler_4.5.0      proxy_0.4-29        withr_3.0.2        
##  [40] doParallel_1.0.17   S7_0.2.1            DBI_1.3.0          
##  [43] RgoogleMaps_1.5.3   ggstats_0.12.0      R.utils_2.13.0     
##  [46] maps_3.4.3          MASS_7.3-65         rappdirs_0.3.4     
##  [49] classInt_0.4-11     gtools_3.9.5        gdsfmt_1.44.1      
##  [52] tools_4.5.0         units_1.0-0         otel_0.2.0         
##  [55] ape_5.8-1           httpuv_1.6.17       R.oo_1.27.1        
##  [58] glue_1.8.0          nlme_3.1-168        promises_1.5.0     
##  [61] grid_4.5.0          cluster_2.1.8.2     reshape2_1.4.5     
##  [64] generics_0.1.4      seqinr_4.2-36       gtable_0.3.6       
##  [67] tzdb_0.5.0          R.methodsS3_1.8.2   class_7.3-23       
##  [70] tidyr_1.3.2         pinfsc50_1.3.0      data.table_1.18.2.1
##  [73] hms_1.1.4           foreach_1.5.2       pillar_1.11.1      
##  [76] stringr_1.6.0       spam_2.11-3         later_1.4.8        
##  [79] splines_4.5.0       tweenr_2.0.3        lattice_0.22-9     
##  [82] PopGenReport_3.1.3  tidyselect_1.2.1    genetics_1.3.8.1.3 
##  [85] knitr_1.51          gridExtra_2.3       xfun_0.61          
##  [88] proto_1.0.0         stringi_1.8.7       yaml_2.3.10        
##  [91] boot_1.3-32         evaluate_1.0.5      codetools_0.2-20   
##  [94] tibble_3.3.1        cli_3.6.5           xtable_1.8-4       
##  [97] Rdpack_2.6.6        jquerylib_0.1.4     dichromat_2.0-0.1  
## [100] Rcpp_1.1.0          pegas_1.4           png_0.1-8          
## [103] StAMPP_1.6.3        parallel_4.5.0      dotCall64_1.2      
## [106] calibrate_1.7.7     gap.datasets_0.0.6  viridisLite_0.4.3  
## [109] mvtnorm_1.3-3       scales_1.4.0        e1071_1.7-17       
## [112] purrr_1.2.1         crayon_1.5.3        combinat_0.0-9     
## [115] rlang_1.2.0

Reading genetic data and clone correction

dat1 <- read.vcfR("md6ct95ms200.vcf")
## Scanning file to determine attributes.
## File attributes:
##   meta lines: 10
##   header_line: 11
##   variant count: 25099
##   column count: 326
## Meta line 10 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 25099
##   Character matrix gt cols: 326
##   skip: 0
##   nrows: 25099
##   row_num: 0
## Processed variant 1000Processed variant 2000Processed variant 3000Processed variant 4000Processed variant 5000Processed variant 6000Processed variant 7000Processed variant 8000Processed variant 9000Processed variant 10000Processed variant 11000Processed variant 12000Processed variant 13000Processed variant 14000Processed variant 15000Processed variant 16000Processed variant 17000Processed variant 18000Processed variant 19000Processed variant 20000Processed variant 21000Processed variant 22000Processed variant 23000Processed variant 24000Processed variant 25000Processed variant: 25099
## All variants processed
dat1
## ***** Object of Class vcfR *****
## 317 samples
## 5870 CHROMs
## 25,099 variants
## Object size: 71.5 Mb
## 0 percent missing data
## *****        *****         *****
## remove clones detected on phylogenetic tree
clones_r <- c(
  "DE1-4","DE1-7","DE1-5",
  "DE2-4",
  "DE3-9","DE3-10","DE3-11","DE3-5",
  "DE4-1","DE4-3","DE4-6","DE4-8",
  "DE5-3","DE5-2","DE5-7","DE5-9",
  "DE7-12","DE7-6","DE7-7","DE7-10","DE7-2","DE7-17","DE7-18","DE7-19","DE7-20","DE7-8",
  "DE9-13","DE9-2r",
  "DE10-11","DE10-12","DE10-2","DE10-5","DE10-7","DE10-8",
  "DE11-12","DE11-15","DE11-14r","DE12-5",
  "DK-4","DK-5",
  "EE-7",
  "FI1-7-1","FI5-2","FI5-10","FI5-8",
  "LT1-9",
  "LT2-3","LT2-4","LT2-5","LT2-6","LT2-7","LT2-1r",
  "LT3-11","LT3-12","LT3-5",
  "LV-12","LV-8",
  "PL2-4","PL2-6-2",
  "PL4-14","PL4-1","PL4-3","PL4-7",
  "PL5-20","PL5-21",
  "PL6-10","PL6-9",
  "PL7-12","PL7-5","PL7-18",
  "PL8-17","PL8-8","PL9-4","PL9-8","PL8-19-",
  "RU1-4","RU2-3","RU2-5","RU2-6",
  "SE3-4-2","SE3-4-2r","SE1-4","SE1-5","SE1-6","SE1-7","SE1-9",
  "UA1-13","UA1-14","UA1-15","UA1-17","UA1-16r",
  "UA2-1","UA2-12","UA2-18","UA2-3")

dat1.pre_gl <- !colnames(dat1@gt)[-1] %in% clones_r
keep <- !colnames(dat1@gt)[-1] %in% clones_r
dat1@gt <- dat1@gt[, c(TRUE, keep)]
dat1.pre_gl <- vcfR2genlight(dat1)
ploidy(dat1.pre_gl) <- 2

dat1.pre_gl <- gl.compliance.check(dat1.pre_gl)
## Starting gl.compliance.check 
##   Processing genlight object with SNP data
##   Warning: data include loci that are scored NA across all individuals.
##   Consider filtering using gl <- gl.filter.allna(gl)
##   The slot loc.all, which stores allele name for each locus, is empty. 
## Creating a dummy variable (A/C) to insert in this slot. 
##   Checking coding of SNPs
##     SNP data scored NA, 0, 1 or 2 confirmed
##   Checking locus metrics and flags
##   Recalculating locus metrics
##   Checking for monomorphic loci
##     Dataset contains monomorphic loci
##   Checking for loci with all missing data
##     Dataset contains loci with all missing dat
##   Checking whether individual names are unique.
##   Checking for individual metrics
##   Warning: Creating a slot for individual metrics
##   Checking for population assignments
##   Population assignments not detected, individuals assigned
##                     to a single population labelled 'pop1'
##   Spelling of coordinates checked and changed if necessary to 
##             lat/lon
## Completed: gl.compliance.check
dat1.pre_gl <- gl.filter.allna(dat1.pre_gl)
## Starting gl.filter.allna 
##   Processing genlight object with SNP data
##   Identifying and removing loci and individuals scored all 
##                 missing (NA)
##   Deleting loci that are scored as all missing (NA)
##   Deleted
##   Deleting individuals that are scored as all missing (NA)
## Completed: gl.filter.allna
dat1.pre_gl@ind.names
##   [1] "DE1-2"   "DE10-1"  "DE10-10" "DE10-3"  "DE11-1"  "DE11-10" "DE11-11"
##   [8] "DE11-14" "DE11-16" "DE11-17" "DE11-2"  "DE11-3"  "DE11-4"  "DE11-5" 
##  [15] "DE11-6"  "DE11-9"  "DE12-1"  "DE12-11" "DE12-12" "DE12-14" "DE12-2" 
##  [22] "DE12-3"  "DE12-4"  "DE12-6"  "DE12-7"  "DE12-8"  "DE12-9"  "DE2-1"  
##  [29] "DE2-2"   "DE2-3"   "DE3-1"   "DE3-2"   "DE3-3"   "DE3-4"   "DE3-8"  
##  [36] "DE4-10"  "DE4-11"  "DE4-2"   "DE4-4"   "DE4-5"   "DE5-1"   "DE5-10" 
##  [43] "DE5-4"   "DE5-6"   "DE6-4"   "DE6-5"   "DE6-6"   "DE6-7"   "DE7-1"  
##  [50] "DE7-11"  "DE7-13"  "DE7-23"  "DE7-25"  "DE7-3"   "DE7-4"   "DE8-1"  
##  [57] "DE8-6"   "DE8-7"   "DE8-8"   "DE8-9"   "DE9-1"   "DE9-10"  "DE9-11" 
##  [64] "DE9-12"  "DE9-15"  "DE9-16"  "DE9-2"   "DE9-4"   "DE9-6"   "DE9-7"  
##  [71] "DE9-9"   "DK-8"    "DK-9"    "EE-5"    "EE-6"    "FI1-1"   "FI1-12" 
##  [78] "FI1-14"  "FI1-15"  "FI1-16"  "FI1-17"  "FI1-4"   "FI1-5-1" "FI1-5-2"
##  [85] "FI1-6"   "FI1-7-2" "FI1-8"   "FI2-10"  "FI2-13"  "FI2-15"  "FI2-17" 
##  [92] "FI2-19"  "FI2-20"  "FI2-3"   "FI2-6"   "FI3-14"  "FI3-17"  "FI3-19" 
##  [99] "FI3-21"  "FI3-24"  "FI3-3"   "FI3-5"   "FI3-9"   "FI4-1"   "FI4-13" 
## [106] "FI4-16"  "FI4-2"   "FI4-22"  "FI4-24"  "FI4-27"  "FI4-3"   "FI4-32" 
## [113] "FI5-1"   "FI5-12"  "FI5-4"   "FI5-5"   "FI5-7"   "FI5-9"   "LT1-10" 
## [120] "LT1-11"  "LT1-4"   "LT1-5"   "LT1-6"   "LT2-1"   "LT2-2"   "LT3-10" 
## [127] "LT3-14"  "LT3-7"   "LT3-9"   "LV-10"   "LV-9"    "PL1-10"  "PL1-5"  
## [134] "PL1-6"   "PL1-7"   "PL1-8"   "PL1-9"   "PL2-10"  "PL2-2"   "PL2-3"  
## [141] "PL2-5"   "PL2-6_1" "PL4-10"  "PL4-11"  "PL4-13"  "PL4-15"  "PL4-2"  
## [148] "PL4-5"   "PL4-8"   "PL5-15"  "PL5-19"  "PL5-4"   "PL6-1"   "PL6-13" 
## [155] "PL6-3"   "PL6-4"   "PL6-5"   "PL6-6"   "PL6-7"   "PL7-10"  "PL7-13" 
## [162] "PL7-14"  "PL7-17"  "PL7-7"   "PL7-8"   "PL7-9"   "PL8-1"   "PL8-11" 
## [169] "PL8-12"  "PL8-14"  "PL8-15"  "PL8-16"  "PL8-18"  "PL8-19"  "PL8-2"  
## [176] "PL8-4"   "PL8-5"   "PL8-6"   "PL8-7"   "PL9-1"   "PL9-10"  "PL9-13" 
## [183] "PL9-5"   "PL9-6"   "PL9-7"   "RU1-1"   "RU1-2"   "RU1-3"   "RU1-5"  
## [190] "RU1-6"   "RU2-1"   "RU2-2"   "RU2-4"   "SE1-1"   "SE1-10"  "SE2-1"  
## [197] "SE2-10"  "SE2-2"   "SE2-3"   "SE2-6"   "SE2-8"   "SE3-1"   "SE3-10" 
## [204] "SE3-2-1" "SE3-2-2" "SE3-4-1" "SE3-5"   "SE3-8-1" "SE3-9"   "UA1-1"  
## [211] "UA1-10"  "UA1-12"  "UA1-16"  "UA1-3"   "UA1-7"   "UA2-11"  "UA2-14" 
## [218] "UA2-16"  "UA2-19"  "UA2-2"   "UA2-5"   "UA2-7"
dat1.pre_gl
##  /// GENLIGHT OBJECT /////////
## 
##  // 222 genotypes,  24,852 binary SNPs, size: 10.4 Mb
##  1012804 (18.36 %) missing data
## 
##  // Basic content
##    @gen: list of 222 SNPbin
##    @ploidy: ploidy of each individual  (range: 2-2)
## 
##  // Optional content
##    @ind.names:  222 individual labels
##    @loc.names:  24852 locus labels
##    @loc.all:  24852 alleles
##    @chromosome: factor storing chromosomes of the SNPs
##    @position: integer storing positions of the SNPs
##    @pop: population of each individual (group size range: 222-222)
##    @other: a list containing: loc.metrics  loc.metrics.flags  verbose  history  ind.metrics
dat1.pre_gl$ind.names
##   [1] "DE1-2"   "DE10-1"  "DE10-10" "DE10-3"  "DE11-1"  "DE11-10" "DE11-11"
##   [8] "DE11-14" "DE11-16" "DE11-17" "DE11-2"  "DE11-3"  "DE11-4"  "DE11-5" 
##  [15] "DE11-6"  "DE11-9"  "DE12-1"  "DE12-11" "DE12-12" "DE12-14" "DE12-2" 
##  [22] "DE12-3"  "DE12-4"  "DE12-6"  "DE12-7"  "DE12-8"  "DE12-9"  "DE2-1"  
##  [29] "DE2-2"   "DE2-3"   "DE3-1"   "DE3-2"   "DE3-3"   "DE3-4"   "DE3-8"  
##  [36] "DE4-10"  "DE4-11"  "DE4-2"   "DE4-4"   "DE4-5"   "DE5-1"   "DE5-10" 
##  [43] "DE5-4"   "DE5-6"   "DE6-4"   "DE6-5"   "DE6-6"   "DE6-7"   "DE7-1"  
##  [50] "DE7-11"  "DE7-13"  "DE7-23"  "DE7-25"  "DE7-3"   "DE7-4"   "DE8-1"  
##  [57] "DE8-6"   "DE8-7"   "DE8-8"   "DE8-9"   "DE9-1"   "DE9-10"  "DE9-11" 
##  [64] "DE9-12"  "DE9-15"  "DE9-16"  "DE9-2"   "DE9-4"   "DE9-6"   "DE9-7"  
##  [71] "DE9-9"   "DK-8"    "DK-9"    "EE-5"    "EE-6"    "FI1-1"   "FI1-12" 
##  [78] "FI1-14"  "FI1-15"  "FI1-16"  "FI1-17"  "FI1-4"   "FI1-5-1" "FI1-5-2"
##  [85] "FI1-6"   "FI1-7-2" "FI1-8"   "FI2-10"  "FI2-13"  "FI2-15"  "FI2-17" 
##  [92] "FI2-19"  "FI2-20"  "FI2-3"   "FI2-6"   "FI3-14"  "FI3-17"  "FI3-19" 
##  [99] "FI3-21"  "FI3-24"  "FI3-3"   "FI3-5"   "FI3-9"   "FI4-1"   "FI4-13" 
## [106] "FI4-16"  "FI4-2"   "FI4-22"  "FI4-24"  "FI4-27"  "FI4-3"   "FI4-32" 
## [113] "FI5-1"   "FI5-12"  "FI5-4"   "FI5-5"   "FI5-7"   "FI5-9"   "LT1-10" 
## [120] "LT1-11"  "LT1-4"   "LT1-5"   "LT1-6"   "LT2-1"   "LT2-2"   "LT3-10" 
## [127] "LT3-14"  "LT3-7"   "LT3-9"   "LV-10"   "LV-9"    "PL1-10"  "PL1-5"  
## [134] "PL1-6"   "PL1-7"   "PL1-8"   "PL1-9"   "PL2-10"  "PL2-2"   "PL2-3"  
## [141] "PL2-5"   "PL2-6_1" "PL4-10"  "PL4-11"  "PL4-13"  "PL4-15"  "PL4-2"  
## [148] "PL4-5"   "PL4-8"   "PL5-15"  "PL5-19"  "PL5-4"   "PL6-1"   "PL6-13" 
## [155] "PL6-3"   "PL6-4"   "PL6-5"   "PL6-6"   "PL6-7"   "PL7-10"  "PL7-13" 
## [162] "PL7-14"  "PL7-17"  "PL7-7"   "PL7-8"   "PL7-9"   "PL8-1"   "PL8-11" 
## [169] "PL8-12"  "PL8-14"  "PL8-15"  "PL8-16"  "PL8-18"  "PL8-19"  "PL8-2"  
## [176] "PL8-4"   "PL8-5"   "PL8-6"   "PL8-7"   "PL9-1"   "PL9-10"  "PL9-13" 
## [183] "PL9-5"   "PL9-6"   "PL9-7"   "RU1-1"   "RU1-2"   "RU1-3"   "RU1-5"  
## [190] "RU1-6"   "RU2-1"   "RU2-2"   "RU2-4"   "SE1-1"   "SE1-10"  "SE2-1"  
## [197] "SE2-10"  "SE2-2"   "SE2-3"   "SE2-6"   "SE2-8"   "SE3-1"   "SE3-10" 
## [204] "SE3-2-1" "SE3-2-2" "SE3-4-1" "SE3-5"   "SE3-8-1" "SE3-9"   "UA1-1"  
## [211] "UA1-10"  "UA1-12"  "UA1-16"  "UA1-3"   "UA1-7"   "UA2-11"  "UA2-14" 
## [218] "UA2-16"  "UA2-19"  "UA2-2"   "UA2-5"   "UA2-7"
length(unique(as.character(dat1.pre_gl@chromosome))) # N of loci
## [1] 5858

Isolation by distance

Preaparing data

p_dist <- bitwise.dist(dat1.pre_gl)
p_dist_mat <- as.matrix(p_dist)

Coordinates <- read_tsv("Coordinates.txt")
## Rows: 311 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): Sample
## dbl (2): Longitude, Latitude
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Coordinates <- as.data.frame(Coordinates)
rownames(Coordinates) <- Coordinates[[1]]   # first column
Coordinates <- Coordinates[,-1] 
Coordinates_no_clones <- Coordinates[ !rownames(Coordinates) %in% clones_r, ]

Coordinates_no_clones <- Coordinates_no_clones[adegenet::indNames(dat1.pre_gl), , drop = FALSE]
stopifnot(identical(rownames(Coordinates_no_clones), adegenet::indNames(dat1.pre_gl)))
identical(rownames(Coordinates_no_clones), adegenet::indNames(dat1.pre_gl))
## [1] TRUE
geo_dist <- distm(Coordinates_no_clones[, c("Longitude", "Latitude")], fun = distHaversine)
max(geo_dist) #maximum geographic distance 
## [1] 1998919
rownames(geo_dist) <- rownames(Coordinates_no_clones)
colnames(geo_dist) <- rownames(Coordinates_no_clones)
shared_samples <- intersect(rownames(p_dist_mat), rownames(geo_dist)) 

gen_dist_matrix <- p_dist_mat[shared_samples, shared_samples]
geo_dist_matrix <- geo_dist[shared_samples, shared_samples]

vegan::mantel(gen_dist_matrix, geo_dist_matrix, method = "pearson", permutations = 10000)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## vegan::mantel(xdis = gen_dist_matrix, ydis = geo_dist_matrix,      method = "pearson", permutations = 10000) 
## 
## Mantel statistic r: 0.1841 
##       Significance: 9.999e-05 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0267 0.0344 0.0410 0.0486 
## Permutation: free
## Number of permutations: 10000

Plot geographic vs genetic distance

gen_vals <- as.vector(gen_dist_matrix[lower.tri(gen_dist_matrix)])
geo_vals <- as.vector(geo_dist_matrix[lower.tri(geo_dist_matrix)])

df <- data.frame(geo = geo_vals/1000, gen = gen_vals)

ggplot(df, aes(x = geo, y = gen)) +
  geom_point(alpha = 0.3, color = "#009E73") +
  geom_smooth(method = "lm", color = "red") +
  labs(x = "Geographic distance (km)", 
       y = "Genetic distance (p-distance)")+
   theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

Isolation by climate

Selection of predictors

bio <- worldclim_global(var = "bio",res = 0.5, path = "climate") # res=0.5=30 arc-sec ≈ 1 km

clim_vals <- terra::extract(bio, Coordinates_no_clones)
anyNA(clim_vals)
## [1] FALSE
clim <- clim_vals[, -1]
colnames(clim) <- paste0("bio", 1:19)
rownames(clim) <- rownames(Coordinates_no_clones)

round(cor(clim), 2) # check correlated parameters
##        bio1  bio2  bio3  bio4  bio5  bio6  bio7  bio8  bio9 bio10 bio11 bio12
## bio1   1.00  0.02  0.60 -0.66  0.44  0.90 -0.72  0.45  0.29  0.70  0.91  0.14
## bio2   0.02  1.00  0.60  0.09  0.45 -0.18  0.32  0.36 -0.47  0.05 -0.08  0.20
## bio3   0.60  0.60  1.00 -0.74  0.12  0.63 -0.56  0.18  0.24  0.06  0.71  0.53
## bio4  -0.66  0.09 -0.74  1.00  0.29 -0.89  0.96  0.12 -0.71  0.06 -0.91 -0.50
## bio5   0.44  0.45  0.12  0.29  1.00  0.06  0.29  0.80 -0.55  0.87  0.08 -0.40
## bio6   0.90 -0.18  0.63 -0.89  0.06  1.00 -0.94  0.14  0.59  0.37  0.99  0.27
## bio7  -0.72  0.32 -0.56  0.96  0.29 -0.94  1.00  0.14 -0.75 -0.06 -0.92 -0.40
## bio8   0.45  0.36  0.18  0.12  0.80  0.14  0.14  1.00 -0.52  0.72  0.17 -0.25
## bio9   0.29 -0.47  0.24 -0.71 -0.55  0.59 -0.75 -0.52  1.00 -0.26  0.56  0.31
## bio10  0.70  0.05  0.06  0.06  0.87  0.37 -0.06  0.72 -0.26  1.00  0.35 -0.34
## bio11  0.91 -0.08  0.71 -0.91  0.08  0.99 -0.92  0.17  0.56  0.35  1.00  0.33
## bio12  0.14  0.20  0.53 -0.50 -0.40  0.27 -0.40 -0.25  0.31 -0.34  0.33  1.00
## bio13  0.09  0.36  0.53 -0.35 -0.28  0.14 -0.23 -0.13  0.12 -0.28  0.21  0.97
## bio14  0.13  0.16  0.52 -0.53 -0.43  0.29 -0.42 -0.26  0.35 -0.37  0.34  0.99
## bio15 -0.22  0.35 -0.31  0.69  0.50 -0.54  0.69  0.30 -0.68  0.31 -0.52 -0.49
## bio16  0.12  0.34  0.54 -0.39 -0.28  0.18 -0.27 -0.14  0.15 -0.27  0.25  0.97
## bio17  0.15  0.14  0.54 -0.56 -0.45  0.31 -0.46 -0.29  0.40 -0.38  0.37  0.99
## bio18  0.18  0.41  0.56 -0.35 -0.14  0.19 -0.23  0.04  0.03 -0.15  0.26  0.93
## bio19  0.04  0.02  0.45 -0.55 -0.60  0.27 -0.47 -0.48  0.52 -0.50  0.31  0.94
##       bio13 bio14 bio15 bio16 bio17 bio18 bio19
## bio1   0.09  0.13 -0.22  0.12  0.15  0.18  0.04
## bio2   0.36  0.16  0.35  0.34  0.14  0.41  0.02
## bio3   0.53  0.52 -0.31  0.54  0.54  0.56  0.45
## bio4  -0.35 -0.53  0.69 -0.39 -0.56 -0.35 -0.55
## bio5  -0.28 -0.43  0.50 -0.28 -0.45 -0.14 -0.60
## bio6   0.14  0.29 -0.54  0.18  0.31  0.19  0.27
## bio7  -0.23 -0.42  0.69 -0.27 -0.46 -0.23 -0.47
## bio8  -0.13 -0.26  0.30 -0.14 -0.29  0.04 -0.48
## bio9   0.12  0.35 -0.68  0.15  0.40  0.03  0.52
## bio10 -0.28 -0.37  0.31 -0.27 -0.38 -0.15 -0.50
## bio11  0.21  0.34 -0.52  0.25  0.37  0.26  0.31
## bio12  0.97  0.99 -0.49  0.97  0.99  0.93  0.94
## bio13  1.00  0.94 -0.28  1.00  0.93  0.98  0.85
## bio14  0.94  1.00 -0.58  0.94  1.00  0.90  0.96
## bio15 -0.28 -0.58  1.00 -0.29 -0.60 -0.24 -0.64
## bio16  1.00  0.94 -0.29  1.00  0.93  0.98  0.85
## bio17  0.93  1.00 -0.60  0.93  1.00  0.88  0.97
## bio18  0.98  0.90 -0.24  0.98  0.88  1.00  0.76
## bio19  0.85  0.96 -0.64  0.85  0.97  0.76  1.00
# I keep: # bio1   – Annual mean temperature
          # bio4   – Temperature seasonality
          # bio12  – Annual precipitation
          # bio15  – Precipitation seasonality

clim_sel <- clim[, c("bio1", "bio4", "bio12", "bio15")]


vif_result <- vifstep(clim_sel, th = 10)
vif_result
## No variable from the 4 input variables has collinearity problem. 
## 
## The linear correlation coefficients ranges between: 
## min correlation ( bio12 ~ bio1 ):  0.1354419 
## max correlation ( bio15 ~ bio4 ):  0.6933374 
## 
## ---------- VIFs of the remained variables -------- 
##   Variables      VIF
## 1      bio1 2.313882
## 2      bio4 4.512395
## 3     bio12 1.484284
## 4     bio15 2.412407
clim_scaled <- scale(clim_sel)
env <- as.data.frame(clim_scaled)

gen_mat <- tab(dat1.pre_gl, NA.method = "mean")
gen_mat <- scale(gen_mat, center = TRUE, scale = FALSE)

Global RDA (redundancy analysis)

rda_env <- rda(gen_mat ~ bio1 + bio4 + bio12 + bio15, data = env)
RsquareAdj(rda_env)
## $r.squared
## [1] 0.09773933
## 
## $adj.r.squared
## [1] 0.0811078
G_rda1<-anova(rda_env, permutations = 10000)
G_rda1
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 10000
## 
## Model: rda(formula = gen_mat ~ bio1 + bio4 + bio12 + bio15, data = env)
##           Df Variance      F    Pr(>F)    
## Model      4   105.67 5.8767 9.999e-05 ***
## Residual 217   975.51                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
G_rda2<-anova(rda_env, by = "term", permutations = 10000)
G_rda2
## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## Model: rda(formula = gen_mat ~ bio1 + bio4 + bio12 + bio15, data = env)
##           Df Variance      F    Pr(>F)    
## bio1       1    37.16 8.2651 9.999e-05 ***
## bio4       1    25.71 5.7185 9.999e-05 ***
## bio12      1    23.74 5.2799 9.999e-05 ***
## bio15      1    19.08 4.2435 9.999e-05 ***
## Residual 217   975.51                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Partial RDA: conditioning on geography

rda_partial <- rda(gen_mat ~ bio1 + bio4 + bio12 + bio15 + Condition(Longitude + Latitude), data = cbind(env, Coordinates_no_clones))
RsquareAdj(rda_partial)
## $r.squared
## [1] 0.05807975
## 
## $adj.r.squared
## [1] 0.04203873
anova(rda_partial, permutations = 10000)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 10000
## 
## Model: rda(formula = gen_mat ~ bio1 + bio4 + bio12 + bio15 + Condition(Longitude + Latitude), data = cbind(env, Coordinates_no_clones))
##           Df Variance      F    Pr(>F)    
## Model      4    62.79 3.5368 9.999e-05 ***
## Residual 215   954.31                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(rda_partial, by = "term", permutations = 10000)
## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
## 
## Model: rda(formula = gen_mat ~ bio1 + bio4 + bio12 + bio15 + Condition(Longitude + Latitude), data = cbind(env, Coordinates_no_clones))
##           Df Variance      F    Pr(>F)    
## bio1       1    29.49 6.6448 9.999e-05 ***
## bio4       1    10.24 2.3071 9.999e-05 ***
## bio12      1    10.91 2.4589 9.999e-05 ***
## bio15      1    12.15 2.7366 9.999e-05 ***
## Residual 215   954.31                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(rda_partial, by = "margin", permutations = 10000)
## Permutation test for rda under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 10000
## 
## Model: rda(formula = gen_mat ~ bio1 + bio4 + bio12 + bio15 + Condition(Longitude + Latitude), data = cbind(env, Coordinates_no_clones))
##           Df Variance      F    Pr(>F)    
## bio1       1    10.87 2.4478 9.999e-05 ***
## bio4       1     8.98 2.0235 9.999e-05 ***
## bio12      1     9.29 2.0921 9.999e-05 ***
## bio15      1    12.15 2.7366 9.999e-05 ***
## Residual 215   954.31                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(rda_partial)$cont$importance[, 1:4]
##                              RDA1        RDA2        RDA3        RDA4
## Eigenvalue            31.12367181 14.69983581 9.676947724 7.294359168
## Proportion Explained   0.03060026  0.01445263 0.009514209 0.007171689
## Cumulative Proportion  0.03060026  0.04505289 0.054567098 0.061738787

Plot partial RDA loadings

var_labels <- c(
  "bio1" = "BIO1 – Annual mean temperature",
  "bio4" = "BIO4 – Temperature seasonality",
  "bio12" = "BIO12 – Annual precipitation",
  "bio15" = "BIO15 – Precipitation seasonality")

site_scores <- scores(rda_partial, choices = 1:2, display = "sites")
site_scores_df <- data.frame(
  Sample = rownames(site_scores),
  RDA1 = site_scores[,1],
  RDA2 = site_scores[,2],
  stringsAsFactors = FALSE)


bp_df <- data.frame(
  Variable = rownames(scores(rda_partial, choices = 1:2, display = "bp")),
  RDA1 = scores(rda_partial, choices = 1:2, display = "bp")[,1] * 4,
  RDA2 = scores(rda_partial, choices = 1:2, display = "bp")[,2] * 4)
bp_df$Label <- var_labels[bp_df$Variable]

eig <- summary(rda_partial)$cont$importance
rda1_var <- round(eig[2,1] * 100, 1)  # proportion explained by RDA1
rda2_var <- round(eig[2,2] * 100, 1)  # proportion explained by RDA2

ggplot() +
  geom_point(data = site_scores_df, aes(x = RDA1, y = RDA2), color = "#009E73", size = 3, alpha = 0.5) +
  geom_segment(data = bp_df, aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
               arrow = arrow(length = unit(0.3,"cm")), color = "black", size = 0.8) +
  geom_text_repel(data = bp_df, aes(x = RDA1, y = RDA2, label = Label), size = 4, fontface = "bold") +
  theme_minimal() +
  labs(
    x = paste0("RDA1 (", rda1_var, "%)"),
    y = paste0("RDA2 (", rda2_var, "%)")) +
  theme(legend.position = "none")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Extract the candidate SNPs

rda_loci <- rda_partial
snp_loadings <- scores(rda_partial, display = "species")
# Environmental variable loadings (biplot)
env_scores <- scores(rda_partial, display = "bp")

library(dplyr)

z_scores <- apply(snp_loadings, 2, function(x) scale(x))  # standardize each axis

# Identify SNPs that are outliers on any axis
outlier_idx <- which(apply(abs(z_scores), 1, max) > 3)

candidate_snps <- snp_loadings[outlier_idx, , drop=FALSE]
loci <- unique(sub("_pos.*", "", rownames(candidate_snps)))

cat("Number of candidate SNPs:", nrow(candidate_snps), "\n")
## Number of candidate SNPs: 1101
cat("Number of candidate loci:", length(loci))
## Number of candidate loci: 671

Plotting RDA1 on the map

# Extract RDA1 and RDA2 scores for each individual
rda_scores <- scores(rda_partial, choices = 1:2, display = "sites")

# Convert to data frame and add sample IDs
rda_scores_df <- data.frame(
  Sample = rownames(rda_scores),
  RDA1 = rda_scores[,1],
  RDA2 = rda_scores[,2],
  stringsAsFactors = FALSE)


# Combine with coordinates (make sure row names match Sample IDs)
plot_data <- cbind(Coordinates_no_clones, rda_scores_df[, -1])
plot_data$Population <- gsub("-.*", "", rownames(plot_data))  # optional: extract population code

# Convert to sf object
plot_data_sf <- st_as_sf(plot_data, coords = c("Longitude", "Latitude"), crs = 4326)
plot_data_proj <- st_transform(plot_data_sf, 3035)

# Extract projected coordinates as plain data frame
plot_data_proj_df <- cbind(
  st_drop_geometry(plot_data_proj),
  st_coordinates(plot_data_proj))

# Reproject both map and data
europe <- ne_countries(continent = "Europe", scale = "medium", returnclass = "sf")
europe_proj <- st_transform(europe, 3035)
plot_data_sf <- st_as_sf(plot_data, coords = c("Longitude", "Latitude"), crs = 4326)
plot_data_proj <- st_transform(plot_data_sf, 3035)

# Extract projected coordinates as plain data frame
plot_data_proj_df <- cbind(st_drop_geometry(plot_data_proj), st_coordinates(plot_data_proj))

ggplot(data = europe_proj) +
  geom_sf(fill = "grey80", color = "grey60", linewidth = 0.8) +
  geom_point(
    data = plot_data_proj_df,
    aes(x = X, y = Y, fill = RDA1),  # use fill for shape 21
    shape = 21,                      
    color = "black",                  # border color
    size = 6,
    stroke = 0.5,                     
    alpha = 1,                      
    position = position_jitter(width = 20000, height = 20000)) +
  coord_sf(xlim = c(3900000, 5700000), ylim = c(2500000, 4500000),
    expand = FALSE) +
  scale_fill_gradient2(
    low = "#009E73", mid ="white", high = "#0072B2",
    midpoint = 0,
    name = "RDA1 score") +
  labs(x = "Longitude", y = "Latitude", title = "Spatial distribution of RDA1 scores") +
  theme_minimal()

scores(rda_partial, display = "bp")
##             RDA1        RDA2
## bio1  -0.6241345 -0.08996384
## bio4   0.1943684 -0.09112245
## bio12  0.6594480 -0.17159814
## bio15 -0.0645662  0.46810385
## attr(,"const")
## [1] 22.10919
climate <-aggregate(clim_sel, by = list(Population = sub("-.*", "", rownames(clim_sel))), FUN = mean)