Testing genomic offset with common gardens in genetically structured black spruce ( Picea mariana)

Black spruce, Picea mariana, genomic offset, Gradient Forest, genomic vulnerability, 16 climate change 17 18 * Corresponding authors: [email protected], nathalie.isabel@nrcan-19 rncan.gc.ca 20 21 22 23 24 25 26 27 28 29 Acknowledgments 30 We thank Marie-Claude Gros-Louis, Esther Pouliot for help with DNA extraction and Jean-31 François Légaré for the preparation of plates. We thank Eric Dussault for collecting samples for 32 DNA sequencing. Climate data access was provided through the Power Analytics and 33 Visualization for Climate Science (PAVICS) platform (Ouranos and CRIM, 2018-2025), which is 34 funded through Ouranos, the Computer Research Institute of Montreal (CRIM), Environment 35 and Climate Change Canada (ECCC), CANARIE, the Fonds Vert and the Fonds d’électrification 36 et de changements climatiques, the Canadian Foundation for Innovation (CFI), and the Fonds 37 de Recherche du Québec (FRQ). This study was funded by the Government of Canada through 38 the Genomics Research and Development Initiative – Genomic Adaption and Resilience to 39 Climate Change (GenARCC) project (2022-2027). 40 41 Conflict of interest disclosure 42 Authors declare no conflicts of interest. 43 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 1

44 Boreal forests play a crucial role in regulating climate via storage and release of carbon. 45 Anticipated changes in climate are predicted to increase mortality and decrease the biomass of 46 many boreal tree species, putting at risk the functioning of this ecosystem and hence its role in 47 carbon absorption and climate change mitigation. Genomic offset methods leverage spatial 48 distribution of genomic diversity and its association with environmental variables to generate 49 maps of populations vulnerable to projected changes in climate. Here we analyse over 60 50 populations and more than 1400 individuals of black spruce ( Picea mariana (Mill.) B.S.P), a 51 dominant boreal forest species, to compare population-level genomic offsets calculated with 52 Gradient Forest with multiple fitness traits measured in four long-term (>40 yr) common 53 gardens. Within common gardens, we found that genomic offset predictions were unaffected by 54 the number or type of markers used for model training, with the exception of LFMM climate-55 associated markers. Most models predicted fitness equally well even when the number of 56 populations in the training set was reduced down to ten or was restricted to a single genetic 57 cluster, suggesting that Gradient Forest can reliably estimate fitness for new populations and 58 novel climates. However, model performances varied among common gardens, with highly 59 accurate fitness predictions in some gardens but contradictory results in others. The accuracy of 60 model predictions was strongly influenced by the choice of climate variables and their 61 relationships with fitness traits. Overall, our results indicate that predicting genomic offsets in 62 genetically structured species across large spatial scales is challenging, because of variation in 63 environmental effects on genotypes and in the interactions among climate variables across the 64 landscape. By capitalizing on our comprehensive validation, we identified the most robust 65 models for projecting black spruce fitness loss in response to future climate change. 66 67 68 69 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 2 70 71 72 73

74 Boreal forests comprise 27% of the total fo rested area on Earth, covering more than 1 billion 75 hectares across the northern hemisphere ( Global Forest Resources Assessment 2020 , n.d.). In 76 Canada, 75% of all forests are located in the boreal zone, making it one of the most important 77 ecosystem and economic service providers, as well as being home to many Indigenous 78 communities in the country. Canadian boreal forests are largely composed of conifer species 79 with very wide distributions, such as black spruce, jack pine ( Pinus banksiana Lamb.), 80 lodgepole pine ( Pinus contorta Dougl.) or white spruce ( Picea glauca (Moench) Voss). These 81 species are well adapted to short-growing seasons and cold temperatures. As such these 82 forests play an important role in the regulation of global climate by storing and releasing carbon. 83 Boreal forests in Canada have been under constant influence of natural disturbances such as 84 extreme climate events, wildfires, insect and disease outbreaks since the retreat of the glaciers 85 about 7000 ya (Brandt et al., 2013). With predicted increase in global temperatures (IPCC, 86 2023), boreal forests are facing intensification of some of those challenges, and increase of 87 drought-related episodes in some parts of its distribution (Price et al., 2013), putting at risk 88 long-term functioning of this ecosystem (Gauthier et al., 2015). It is becoming increasingly 89 urgent to understand how boreal species will respond to these challenges, and test what are the 90 potential benefits of assisted migration to species and tree populations at risk of maladaptation. 91 Ecological methods used to predict forest species dynamics in changing environments either 92 model species distributions by examining shifts in the suitability of their habitats across the 93 landscape (Prasad et al., 2020) or model biomass losses by examining species-specific 94 functional traits such as tree mortality caused by drought, or migration failure (Aubin et al., 2018; 95 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 3 Liu et al., 2023). These methods do not take into account the evolutionary processes shaping 96 intraspecific genetic variation across the landscape, which can lead to divergent populations’ 97 responses to the same climate variables. Standing genetic diversity within populations has been 98 shaped by past demographic events, such as local extinctions or expansions, by local selective 99 pressures as well as by interactions with other species. These factors ultimately determine the 100 ability of populations to re spond to changing environments and can enhance the predictive 101 accuracy of species distribution models (Aguirre-Liguori et al., 2021; Urban et al., 2016). Conifer 102 tree species often display local adaptation to environmental conditions, as demonstrated with 103 common garden experiments, where population-level functional traits often shift when growing 104 conditions differ from those found at the seed source (Depardieu et al., 2020; Guo et al., 2022; 105 Housset et al., 2018; Lortie & Hierro, 2022; Oleksyn et al., 1998). Genomic studies further 106 support this, revealing allele frequency changes in genes associated with climate-related traits 107 (Bergmann, 1978; Depardieu et al., 2021; Evans et al., 2014). An increasing number of studies 108 takes into account population-level variation in phenotypes or genotypes in order to model 109 species productivity predictions (Bradley St Clair & Howe, 2007; Isaac-Renton et al., 2014; 110 Patsiou et al., 2020; Robert et al., 2024). 111 Another group of predictive methods aims to estimate genotype–environment mismatch 112 (genomic offset) under climate change using species-wide genomic variation. Approaches such 113 as Gradient Forest (Ellis et al., 2012; Fitzpatrick & Keller, 2015), Risk of Non-Adaptedness 114 (Rellstab et al., 2016), Latent Factor Mixed Models (Gain & François, 2021), or Redundancy 115 Analysis (Capblancq & Forester, 2021) model linear or non-linear relationships between genetic 116 variants and environmental variables across the landscape. These models are then used to 117 predict genomic offset under future climate scenarios. Genomic offset is typically calculated as 118 the Euclidean distance between the genomic composition of a genotype under current 119 conditions and its projected composition under new conditions, and it has been theoretically 120 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 4 linked to changes in fitness when a polygenic trait is under stabilizing selection toward an 121 optimal phenotype (Capblancq et al., 2025; Gain et al., 2023). Such methods come with their 122 own limitations; for instance the underlying assumption is that genotypes are most fit in their 123 local climate, and that maladaptation to new conditions is instantaneous. Therefore, even 124 though these methods hold promise for straightforward and phenotype-free risk assessment of 125 species to changing climate (Rellstab et al., 2021), validation studies are needed to assess their 126 utility for conservation or management purposes in different contexts. 127 To validate genomic offset methods, several studies have conducted either simulations (Gain et 128 al., 2023; Láruson et al., 2022; Lind & Lotterhos, 2025a, 2025b) or empirical validation in 129 common gardens in which anticipated climate change in some future period is substituted by 130 different climate conditions experienced by the genotype transplanted to a common garden 131 experiment (Archambeau et al., 2025; Capblancq & Forester, 2021; Fitzpatrick et al., 2021; Gain 132 et al., 2023; Lachmuth et al., 2024; Lind et al., 2024; Rhoné et al., 2020; Verrico et al., 2025). 133 The overall consensus across these studies is that genomic offset is often correlated with 134 fitness, and often predicts fitness better than climate or geographical distance alone, however 135 different models can produce inconsistent patterns of maladaptation (Lind et al., 2024) and 136 correlations vary depending on the studied species, data- and garden-specific factors 137 (Fitzpatrick et al., 2025). Simulations demonstrated that accuracy of genomic offset is higher 138 when sampling is more dense along environmental gradients (Láruson et al., 2022), and when 139 selection is strong (Lind & Lotterhos, 2025b). On the other hand, differences in effective 140 population sizes (Láruson et al., 2022), presence of novel environments, or presence of non-141 causal environmental variables can decrease accuracy of genomic offset predictions (Lind & 142 Lotterhos, 2025b). Finally, the choice of markers, whether they were random, neutral, or 143 adaptive, either had limited impact on genomic offset estimates (Fitzpatrick et al., 2021; Láruson 144 et al., 2022; Lind et al., 2024; Lind & Lotterhos, 2025b) or improved the accuracy in other 145 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 5 studies (Gain et al., 2023). These findings highlight the importance of exploring the limitations 146 of genomic offset across different species, with particular focus on unsolved questions, such as 147 understanding the impact of genetic structure on the variability of genomic offset patterns, 148 determining which fitness traits are best representing local adaptation, and understanding what 149 impact has the choice of the environmental variables, or of the sampling scheme, on the 150 accuracy of the offset predictions. 151 Our study focuses on black spruce, a widespread boreal species extending from Alaska to 152 Newfoundland and Labrador (Burns & Honkala, 1990). This dominant component of Canadian 153 boreal forests displays adaptation to variable soil and climate conditions, ranging from poorly 154 drained bogs to the northern treeline in northwestern Canada, where it can be found in cold, dry 155 climates and on permafrost-underlain soils. Historically, black spruce has reached its current 156 range from at least two to three distinct glacial refugia situated in the south after the retreat of 157 the glacier (Gérardi et al., 2010; Jaramillo-Co rrea et al., 2004). Studies using mitochondrial, 158 cytoplasmic and nuclear markers identified three main genetic clusters, roughly covering the 159 western, central and eastern part of Canada (Gérardi et al., 2010; Girardin et al., 2021; Prunier 160 et al., 2012), suggesting that current distribution was formed via multiple expansions from the 161 southern refugia, and indicating that black spruce range is composed of partially independent 162 genetic pools in which climate adaptation occurs (Prunier et al., 2012). 163 The impact of climate on the growth of black spruce has been extensively studied, in particular 164 using dendroecological data. Interplay between temperature and precipitation determines water 165 availability stress which is one of the most important factors limiting growth, and therefore 166 productivity across the black spruce range (Girardin et al., 2024). Black spruce also displays 167 generally poor recovery after drought stress compared to other species, such as jack pine 168 (Marchand et al., 2021, 2025). The western and central parts of black spruce distribution are 169 characterised by climate where the ratio of available moisture supplied by precipitation to 170 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 6 evapotranspiration demands is much lower than in the east, limiting tree growth (Hogg et al., 171 2013; Price et al., 2013). While higher temperatures can increase the growth rate of some 172 populations in the northwest, often the opposite is true, and temperature-induced droughts have 173 negative effects on growth (Sniderhan et al., 2021; Walker et al., 2015). Growth of eastern more 174 humid black spruce populations is on the other hand mostly limited by cold temperatures in the 175 north, however, higher temperatures can induce water availability constraints as well as higher 176 respiration costs in the south (Girardin et al., 2016). 177 Taken together, forecasts very broadly predict that increasing temperatures and droughts will 178 lead to higher vulnerability espec ially of western and central North American black spruce 179 populations where precipitation is the limiting factor (Lesven et al., 2024). While northeastern 180 populations can potentially experience growth increase in the near future, this effect is expected 181 to be transient (D’Orangeville et al., 2018). So far only few studies explicitly took into account 182 genetic variation of black spruce, focusing on modelling productivity and resistance to extreme 183 drought events in the context of assisted migration both on the latitudinal and longitudinal scale 184 (Girardin et al., 2021; Robert et al., 2024, 2025). These studies showed that genetic divergence 185 between black spruce genetic clusters has an effect on phenotypic response to climate that can 186 accumulate over time (Girardin et al., 2021; Robert et al., 2024). Models suggested that 187 northern-most sites may experience transient increase in growth due to decrease in cold 188 limitation, however eventually all populations were predicted to decrease in biomass with the 189 western cluster exhibiting the lowest growth compared to the eastern and central clusters 190 (Robert et al., 2024). It is important to note though that other factors linked to climate can 191 exacerbate or offset these predictions. This includes outbreaks of insects such as eastern 192 spruce budworm (Bellemin-Noël et al., 2021), decreas ed ability to recover from fires (Baltzer et 193 al., 2021), loss of protective snow cover (Marquis et al., 2022), changes in tree species 194 composition, which in turn can influence competition as well as composition of soil symbiotic 195 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 7 partners, and changes in bud phenology, where higher temperatures can lead to earlier bud 196 break in the spring and higher potential frost damage (Guo et al., 2022; Mura et al., 2022; Rossi, 197 2015). While in this study we focus specifically on direct climate impact, the common garden 198 experiments of black spruce provide ideal conditions to test how genome-wide variants can add 199 nuanced resolution to predicting tree response to changing climate. 200 Here we examined >1,400 black spruce trees for genotypes and 17 phenotypic traits originating 201 from 66 populations and grown across four common gardens. We leveraged nearly 30,000 202 DArTseq single nucleotide variants, in combination with climate and phenotypic information to 203 validate Gradient Forest method, and predict genomic offset of black spruce to future climate 204 across the species distribution. Our main objectives were: i) to determine fine-scale genetic 205 structure and distribution of genetic diversity across the black spruce range, ii) to validate 206 Gradient Forest genomic offset using a space-for-time substitution based on the phenotypes 207 produced in the common garden experiment and iii) to predict genomic offset of black spruce to 208 projected climate. Our extensive validation of genomic offset in a species with a pronounced 209 genetic structure provides valuable insights for informing future forestry management strategies 210 in the face of changing climate. 211

212 Sample Collection and DNA Extraction 213 Samples were collected between 2014 and 2019, from 70 provenances of black spruce, to 214 which we refer to as populations throughout the text. A total of 2,237 specimens were collected, 215 encompassing the species' entire natural range in Canada and the USA (Fig. 1A, Suppl. Table 216 1, Suppl. Dataset 1). These samples were gathered from common gardens located in 217 Chibougamau (CH, Quebec), Acadia (AC, New Brunswick), Mont-Laurier (ML, Quebec), Peace 218 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 8 River (PR, Alberta), and Valcartier (VL, Quebec), all of which were initially established in 1974 219 or 1975 as part of the Range-Wide Provenance Study, with seeds sown in 1970 ((Morgenstern, 220 1978); Suppl. Table 2). Only a small fraction of samples (n=37, nine of which remained after 221 filtering) was collected from VL and used exclusively for genotyping. Additionally, we collected 222 samples from 37 red spruce ( P. rubens Sarg.) trees, from various populations within the 223 allopatric zone of the species (Suppl. Dataset 2). Red spruce is known to hybridize with black 224 spruce, we therefore included it to account for genetic admixture between the two species. 225 The DNA required for sequencing was extracted from 30-50 mg of either frozen needle tissue 226 (for the CH, ML and VL sites) or cambial tissue (for AC and PR sites) with a Nucleospin 96 Plant 227 II kit (Macherey-Nagel, Bethlehem, PA) using the centrifugation processing protocol with a cell 228 lysis step with PL2 buffer for 1h at 65°C. Cambium (with phelloderm) tissue was obtained from 229 bark samples collected with a 1 cm diameter punch sterilized with 70% alcohol between each 230 sample. Prior to extraction, all tissues were ground to powder with a Mixer Mill MM300 (Retsch 231 GmbH, Haan, Germany) after being plunged into liquid nitrogen for two minutes (this process 232 was repeated twice). 233 Sequencing and Genotyping 234 For single nucleotide polymorphisms (SNPs) discovery, we employed the DArTseq™ method, 235 which shares similarities with genotyping-by-sequencing (GBS) but includes a complexity 236 reduction step targeting low-copy sequences within the genome. Our procedure involved 237 digesting genomic DNA with PstI and MseI restriction enzymes, ligating barcoded adapters, 238 amplifying the resulting products via PCR, and sequencing them using a HiSeq 2500 system. 239 This entire process was conducted by Diversity Arrays Technology, headquartered in Bruce, 240 Australia, as described by (Kilian et al., 2012). Ou r SNP analysis focused exclusively on biallelic 241 SNPs for subsequent investigation. The DNA plates containing extracted samples were sent 242 and processed in three consecutive years: 2018, 2019, and 2020, referred to herein as batches 243 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 9 1, 2, and 3 (Suppl. Dataset 1). Genotypes were generated separately for each of these three 244 batches, and at the conclusion of the project, all sequences produced were merged and 245 clustered together to create the final dataset for analysis. 246 Data processing 247 Genetic data 248 Single nucleotide variants (SNPs) were generated from ~70 nt sequence clusters by setting 249 different sequence identity thresholds producing four independent sets of variants. After 250 checking SNP genotyping rates and ensuring no overlap among sequence tags from the four 251 files, SNP sets resulting from the three lowest clustering thresholds were combined together. 252 SNPs were filtered separately for the black spruce and red spruce samples. We filtered out 253 variants with over 50% of missing genotypes, heterozygosity above 50% and the minor allele 254 frequency (MAF) below 1%. We run admixture analysis to determine population structure and at 255 the same time identify potential mislabelled samples. The number of genetic clusters were 256 determined with ADMIXTURE v1.3.0 (Alexander et al., 2009) using all black and red spruce 257 individuals. First genotypes were saved to plink format using the function gl2plink() from dartR 258 package v2.9.7 (Mijangos et al., 2022). Genetic ancestries were assigned for the K number of 259 clusters, where K ranged between 2 and 15. 5-fold cross-validation (cv) error was estimated 260 with each run. CV error curve flattened starting with K=4, pointing to three main black spruce 261 genetic clusters (West, Central, East) and one red spruce cluster. Several samples showed 262 assignment to groups from distant clusters, implying potential mislabelling of samples. For the 263 final dataset we removed the most dubious samples, including 5 samples in the region A and B 264 (Alaska) assigned to the Central or East genetic cluster, and 1 sample from the region I 265 (Newfoundland) assigned to the West cluster. In addition we removed samples with a call rate 266 below 85% retaining 1467 black spruce samples from 66 populations (mean of 22 samples per 267 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 10 population, range 1-49). A dataset of uniquely black spruce variants consisted of 29,978 biallelic 268 SNPs. 269 270 271 Climate data 272 We used the BioSim v11 (Fortin et al., 2022) to simulate daily weather at the populations’ 273 origins and that of the common gardens for the period 1961-1990, then we calculated 43 climate 274 variables representing temperature, precipitation and the interplay between the latter such as 275 climate moisture index (CMI) at different temporal scales (at the annual or seasonal scale; see 276 Suppl. Table 3, Suppl. Dataset 3 for a complete list of variables and the climate aggregate over 277 the full period for 70 populations and 5 common gardens). We obtained simulated daily climate 278 data between 2021 and 2100 for the three emission scenarios (shared socio-economic 279 pathways): SSP2-4.5, SSP3-7.0 and SSP5-8.5, using an ensemble of 13 climate models from 280 the Ouranos ESPO-G6-R2v1.0.0 dataset (Lav oie et al., 2024). Annual and seasonal climate 281 variables corresponding to a subset of those obtained for the past period (a total of 23) were 282 calculated from daily climate data. Mean annual CMI for the future period was obtained by 283 importing daily variables to BioSim. To ensure that variables obtained with BioSim and those 284 estimated from PAVICS datasets were comparabl e, we visually inspected the temporal 285 distributions of annual climate means per population to confirm that means calculated for past 286 and future periods followed a continuous trend and retained 14 of them (CMI, mayMinT, TP, 287 MAT, MMinT, MMaxT, sumTP, WSTP, winMT, sprMT, sumMT, fallMT, MCMT, MWMT). For 288 clarification we refer to these 14 climate variables as a “small climate set”. 289 Population structure analyses 290 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 11 Population structure was determined using an ADMIXTURE program as described above and 291 principal component analysis was used to visualize ordination of genetic clusters. Genotypes 292 with the black spruce individuals were saved to gds format using the function gl2gds() from 293 dartR package v2.9.7 (Mijangos et al., 2022). PCA was run with snpgdsPCA() from R package 294 SNPRelate v1.40.0 (Zheng et al., 2012). 295 To investigate spatial structure we considered 63 black spruce populations with at least five 296 individuals. Weir and Cockerham estimate of FST between pairs of populations was calculated 297 with fst_WC84() function in R package assigner v0.6.0 (Gosselin et al., 2019). To do that, 298 genotypes were first transformed to the tidy object using R package radiator v1.3.5 (Gosselin, 299 2020). To obtain a matrix of distances in km between populations, coordinates were first 300 transformed to a spatial object sf using a WGS 84 coordinate reference system with st_as_sf() 301 function from R package sf v1.0-19 (Pebesma, 2018). Distances were calculated with the 302 function st_distance() from R package sf v1.0-19 (Pebesma, 2018). To test if population genetic 303 divergence increased with distance, a Mantel test was calculated by correlating the matrices of 304 genetic distances ( FST - (1 - F ST)) and geographic distances with the function mantel() from R 305 package vegan v2.6-8 (Oksanen et al., 2025), with 9999 permutations, and the spearman rank 306 correlation method. In parallel, to test for spatial autocorrelation and correlation between genetic 307 and spatial distribution, we generated Moran’s eigenvector maps (MEMs). First a spatial 308 weighting matrix was created for populations with function chooseCN() and parameters “type = 309 5, d1 = 0, d2 = 20” from adegenet v2.1.10 (Jombart, 2008) and nb2listw() from spdep v1.3-8 310 (Bivand & Wong, 2018). MEMs were obtained with function mem() from package adespatial 311 v0.3-24 (Guénard & Legendre, 2022). To conduct spatial PCA, we estimated allele frequencies 312 within populations with function gl.alf() from package dartR v2.9.7 (Mijangos et al., 2022) and 313 removed SNPs with missing frequencies in any of the populations. sPCA was run with the 314 spca() function from the adegenet package v2.1.10 (Jombart, 2008), using a spatial weighting 315 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 12 matrix estimated above. Moran’s coefficient was calculated with the moran.randtest() function 316 from package adespatial v0.3-24 (Guénard & Legendre, 2022) to test if genetic variation 317 represented with principal components is spatially autocorrelated. 318 Pairwise FSTs between clusters were estimated in the same way as between the populations. 319 Individuals were assigned to one of the five clusters identified with admixture with the number of 320 clusters K=6 excluding red spruce (West, Central, East, WI and ME). Each individual was 321 assigned to a cluster based on the dominant ancestry assignment. Analysis of molecular 322 variance (AMOVA) was conducted at three levels: cluster, population and individual. Only 323 populations with at least 5 individuals were retained. AMOVA was estimated with the function 324 poppr.amova() from R package poppr v2.9.6 (Kamvar et al., 2014). Heterozygosity was 325 calculated individually for each population, using the original dataset only filtered for variants not 326 sequenced in more than 50% of samples, but with rare alleles retained. We focused on 60 327 populations with seven or more trees, and randomly drew seven genotypes five times. For each 328 subset, we removed any missing and monomorphic sites, and calculated observed 329 heterozygosity as the number of heterozygous SNPs for each individual divided by all sites, and 330 averaged across samples in each population. Expected heterozygosity was calculated for each 331 population as the sum of squared frequency of allele 1 and 2 in each site, averaged across sites 332 and subtracted from 1. 333 Spatial associations between genotypes and climate 334 To visualize clustering of samples by climate, we run PCA on 43 scaled and centered climate 335 variables using prcomp() function from the base R package stats. To estimate how much 336 distribution of genetic diversity is correlated with climate we conducted partial redundancy 337 analysis (pRDA) by fitting three main climate principal components (PCs), three main Moran 338 eigenvectors (MEMs), and two main PCs from genetic diversity ordination to a matrix of 339 population allele frequencies. Models were fitted with the function rda() to determine the impact 340 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 13 of genetic, spatial and climate explanatory variables, one class (spatial, genetic, climate) at a 341 time while controlling for the remaining variables. Model R2 was estimated with the function 342 RsquareAdj() and the significance of correlation of each class of variables with genetic variation 343 was tested with a permutation test implemented in the function anova.cca(). All functions are 344 found in the vegan package v2.6-8 (Oksanen et al., 2025). 345 To determine which climate variables are best associated with genetic variation we run Gradient 346 Forest (Ellis et al., 2012; Fitzpatrick & Keller, 2015) using 43 climate variables and 3 MEMs 347 depicting spatial patterns of distribution. First, we used the allele frequency matrix generated in 348 previous steps for 63 populations, removed any non-variant sites and selected a random subset 349 of 1000 sites. We fitted the model using the gradientForest package v0.1-37 (Ellis et al., 2012), 350 with 500 bootstrapped trees and correlation threshold of 0.5. Variables were ranked based on the 351 mean raw accuracy importance. 352 Testing genomic offset with common gardens 353 We used common garden experiments to test the accuracy of genomic offset predictions 354 calculated with Gradient Forest. We chose this method because previous studies showed that it 355 performs equally well or better compared to other methods (Capblancq & Forester, 2021; Lind & 356 Lotterhos, 2025b). To validate Gradient Forest, in each common garden, we trained the model 357 using population-specific climate and allele frequencies, and then calculated population 358 genomic offsets projected to common garden climate conditions. Next we conducted 359 Spearman’s rank correlation between genomic offsets, proxies of population maladaptation to 360 common garden conditions, and population phenotypic trait means, proxies for population 361 fitness, expecting more negative correlation coefficients for better-performing models. To train 362 the models, individual trees were split in each population into test and train datasets, to bypass 363 the issue of testing model accuracy for the same genotypes on which the model was trained. To 364 test which factors influence model performance, several aspects of the data were explored, 365 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 14 including type and number of markers, scope of the training datasets (the number of 366 populations, and the number of genotypes per population, the genetic ancestry of the training 367 dataset), scope of the test datasets, and the choice of climate parameters. Details of each step 368 are described below. 369 370 Phenotypic data 371 Population means and standard errors of all traits used as a proxy for fitness are given in Suppl. 372 Dataset 4. Biomass increments and wood densities (PR: n=397, ML: n=612, CH: n=507, AC: 373 n=481) were obtained from wood cores described and published in (Robert et al., 2024). Tree 374 biomass increment (kg) was calculated as the sum of increments between years 1980 and 375 2015, and in eight non-overlapping five-year intervals between 1980 and 2020. Wood density 376 was reported as the average wood density for the ring from the year 2015. Height in cm and 377 diameter at breast height (DBH) in cm, were measured in 2019 (PR: n=400), 2015 (ML: n=667), 378 2016 (CH: n=451), and 2017 (AC: n=492). Survival was reported as the % of surviving trees per 379 population in 2011 (PR: n=26), 2007 (ML: n=42), 2006 (CH: n=38), and 2003 (AC: n=40). Four 380 indices of resilience to extreme drought events (sensu (Lloret et al., 2011)) were calculated from 381 the wood cores and averaged over 14 extreme drought events occurring in one of the four 382 gardens (Robert et al. in prep). These indices included resistance (Rt) - growth during the year 383 of the drought over the average growth of the two years prior to the drought (higher ratio higher 384 resistance), recovery (Rc) - the average growth of the two years after the drought over the 385 growth during the year of the drought (higher ratio higher recovery), resilience (Rs) - the ratio of 386 growth after the drought event over the growth before the event (higher ratio higher resilience), 387 and relative resilience (Rr) - resilience minus resistance (higher ratio higher relative resilience). 388 To control for random effects in the common gardens, we applied best linear unbiased 389 prediction (BLUP) to each trait except for survival, using the following model: 390 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 15 /g1861/g1866/g1856/g1861/g1874/g1861/g1856/g1873/g1853/g1864 /g1872/g1870/g1857/g1857 /g1872/g1870/g1853/g1861/g1872 ~ /g1854/g1864/g1867/g1855/g1863 /g3397 /g46661|/g1868/g1867/g1868/g1873/g1864/g1853/g1872/g1861/g1867/g1866/g4667 /g3397 /g46661|/g1854/g1864/g1867/g1855/g1863: /g1868/g1867/g1868/g1873/g1864/g1853/g1872/g1861/g1867/g1866/g4667. Biomass increments were log transformed. Blocks with less than 10 individuals were excluded. 391 Corrected traits were then averaged for each population. 392 393 Climate transfer distance 394 To investigate local adaptation to climate, a Spearman’s rank correlation was calculated 395 between phenotypic traits and climate transfer distance while accounting for spatial 396 autocorrelation using modified t test with a function modified.ttest() from R package SpatialPack 397 v0.4-1 (Vallejos et al., 2020). P-values were corrected for multiple testing with the Benjamini-398 Hochberg procedure at a false discovery rate of 0.05. Climate transfer distance was calculated 399 as the Euclidean distance of each population to its origin in climate space defined by 43 climate 400 variables. To calculate correlations, populations with less than four trait values were removed. 401 Defining Training and Test Sets 402 To compare population genomic offsets with phenotypic traits, different samples were selected 403 for training the model (train dataset) and scoring phenotypic traits (test dataset). In model 404 validation, we ignored two populations containing samples classified into ME or WI genetic 405 clusters, to minimize the potential confounding effects from genetically distinct locations. An 406 automated script was implemented to partition individuals into train and test datasets. First, all 407 samples, which had only phenotypic information (low quality genotypes were filtered out), were 408 put in the test dataset, since these could not be used for training Gradient Forest. Next, the rest 409 of the samples were put in the train datasets. If a population was absent from the test dataset or 410 comprised less than three samples, samples were transferred from the corresponding 411 population in the train dataset, unless the train population comprised less than five samples. In 412 this way we aimed to have at least three samples for estimating mean phenotypic traits per 413 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 16 population in the test dataset, and at least five genotypes per population for training, although 414 occasionally smaller sample sizes could occur. In the case of AC, many genotypes were 415 missing due to low quality, therefore test datasets had more populations, and more samples 416 than train datasets. Selection of train and test datasets and following analyses were done 417 separately for each trait, as the number of samples and populations varied across them. This 418 accounted for 64 datasets (16 traits excluding survival and four gardens). Numbers of 419 populations and trees per garden and trait selected for test and train datasets are given in the 420 Suppl. Tables 4-5. 421 Detecting Climate-Associated Variants 422 To detect genetic variants showing associations with distinct climate niches (precisely with the 423 three main PC axes on 43 climate variables) we used three approaches: LFMM, RDA, and RDA 424 with correction for population structure. Genotype-environment associations for each SNP were 425 tested using train datasets combined across gardens. LFMM v2 was run with lfmm2() from 426 package LEA v3.18 (Gain & François, 2021), using K=3 latent factors (corresponding to three 427 black spruce genetic clusters). P-values were calculated with lfmm2.test() to test the 428 significance of associations between each SNP and a combined effect of three climate PCs. P -429 values were then adjusted with Bonferroni multiple test correction. RDA was conducted with 430 rda() from package vegan v2.6-8 (Oksanen et al., 2025), first with the three climate PCs as only 431 explanatory variables, and second by additionally conditioning for population structure using the 432 two main axes from genetic variation ordination. Variant P-values were obtained by calculating 433 Mahalanobis distances of each variant from the centroid of RDA variant loadings and comparing 434 them with chi-squared distribution using the method recipe provided in (Capblancq et al., 2018), 435 and finally adjusting with Bonferroni multiple test correction. To estimate the amount of genetic 436 variation represented by the selected list of outliers, observed heterozygosity was calculated. 437 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 17 Heterozygosities per population were estimated on a combined set of all available genotypes 438 from the four gardens and without filtering for missing data. 439 Modeling with Gradient Forest 440 We tested the performance of genomic offset predictions using models trained on each garden 441 train dataset. Genotype data was first prepared by retaining only unique variants per sequence 442 tag, calculating allele frequency per population separately for each common garden, and 443 removing any invariant sites or sites with missing data. We tested several options of input 444 genotypes, including a random subset of 100 (0.1K), 1000 (1K) or all variants (~29K), subsets of 445 variants not filtered for minor allele frequency (01.K-lf, 1K-lf) as well as variants associated with 446 climate (outliers detected with LFMM, RDA and RDA-struct when controlling for population 447 structure). Gradient Forest was run using R package gradientForest v0-1.37 (Ellis et al., 2012), 448 with 500 trees and a correlation threshold 0.5. Genomic offset was calculated as the Euclidean 449 distance between vectors of allelic turnover functions projected for populations in their native 450 climate and climate specific to common gardens where they were growing. Climate was defined 451 by the three main PC axes obtained from 43 climate variables. 452 Sensitivity analysis 453 We explored the impact of downsampling of the train or test datasets, and the impact of the 454 selection of climate variables on genomic offset predictions. For this purpose, all models were 455 trained using 1000 random variants. To test the impact of downsampling of test and train 456 datasets, we systematically removed i) individual populations from the train dataset, leaving 457 between 6 to 20 populations (with a step of 2) randomly selected in five draws, ii) individual 458 genotypes from populations from the train dataset, leaving between a maximum of 2 to 10 459 genotypes per population (with a step of 2) random ly selected in ten draws, and iii) individual 460 genotypes from populations from the test dataset, leaving between a maximum of 2 to 10 461 genotypes per population (with a step of 2) randomly selected in ten draws. In addition, we 462 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 18 tested model performance when models were trained on populations coming from a single 463 genetic cluster (West, Central or East). To test the impact of the choice of climate variables on 464 model performance, we trained models using two sets of climate variables: i) the three main 465 principal component (PC) axes derived from 43 climate variables, and ii) top five variables (dPP, 466 CMI, TP, WSTP, and fallMT) identified by Gradient Forest as most important from a smaller 467 climate set supplemented with a photoperiod variable (n=15). 468 To validate genomic offsets, Spearman’s rank correlation was calculated if at least five 469 populations had estimates of both genetic offset and mean trait value. We also ensured that at 470 least five trees per population were present in the train dataset and three trees per population 471 were present in the test dataset. 95% confidence intervals were computed with 1000 bootstraps. 472 To calculate correlation coefficient within clusters, each population was classified into one of the 473 three clusters based on individual-ancestry frequency information. First, trees were assigned to 474 a genetic cluster, and then populations. Because rho correlation coefficients between fitness 475 and climate transfer distance were very similar regardless of whether they were calculated with 476 or without taking into account spatial autocorrelation, therefore for consistency, in the validation 477 part we report all Spearman rho estimates estimated without correction. To investigate 478 consistency in genomic offset predictions between models we calculated Pearson’s correlation 479 coefficients. 480 Genomic vulnerability of black spruce to future climate 481 Using a small climate set (CMI, mayMinT, TP, MAT, MMinT, MMaxT, sumTP, WSTP, winMT, 482 sprMT, sumMT, fallMT, MCMT, MWMT), future simulated annual means year were averaged 483 across four 20-year periods: 2022-2040, 2041-2060, 2061-2080, and 2081-2100. Current (1960-484 1990) and future-period climate variables for three emission scenarios were combined together, 485 scaled, and clustered with KMeans() function in the scikit-learn v1.7.0 python package, using 486 n_init=’auto’ parameter and the number of clusters ranging from 1 to 25. We then implemented 487 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 19 the elbow method with the KElbo wVisualizer() function from the yellowbrick v1.5 python 488 package, to select the optimal number of clusters. The method measures when adding another 489 cluster does not considerably decrease the distortion score (measure of the tightness of data 490 points to their cluster). The method split all data points into seven clusters, which were 491 visualised on the two main principal components obtained for 14 climate variables using PCA() 492 function in the scikit-learn v1.7.0 python package. 493 Gradient Forest was run as described above on a subset of 1000 randomly selected SNPs and 494 allele frequencies calculated for each population including those from clusters ME and WI to 495 span the full range of black spruce genetic diversity. SNPs with any missing frequency 496 estimates or no variation were excluded. Genomic offsets were calculated for the four 20-year 497 periods and three emission scenarios mentioned above. Two models were trained: i) model 498 trained on all populations and the three main principal components derived from a small climate 499 set (14 variables mentioned above and photoperiod) to project genomic offsets in the 500 Central/East cluster, and ii) model trained on populations from the West cluster and a subset of 501 five selected variables (dPP, CMI, TP, WSTP, and fallMT), to project genomic offsets in the 502 West cluster. 503 All analyses were performed with R v4.4.2 (R Core Team, 2021) and python v3.12. 504

505 Population structure of the black spruce 506 We analyzed 1467 black spruce individuals from 66 populations (Fig. 1A, Suppl. Table 1, Suppl, 507 Dataset 1). We identified 3 main genetic clusters of black spruce, West, Central and East, with 508 most samples showing mixed ancestry between West and Central or Central and East clusters 509 (Fig. 1B-C) deriving from either gene flow or incomplete lineage sorting. PCA showed similar 510 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 20 results, with the 1st PC separating samples from cluster West, and the 2nd PC spreading 511 samples from clusters East and Central, but with no clear separation (Fig. 1D). With the higher 512 K number of clusters, the admixture program distinguished two additional subclusters: WI 513 (samples from a population in Wisconsin) within the Central cluster, and ME (samples from a 514 population in Maine) within the East cluster (Fig. 1B-C). FST was 0.009-0.012 between the 515 cluster West and the clusters East or Central, about ten times higher than between clusters 516 Central and East (FST=0.001, Fig. 1E). Clusters WI and ME were most differentiated from cluster 517 West ( FST between 0.062 and 0.067) and also from each other ( FST=0.032). AMOVA showed 518 low but significant variance between clusters (% variance between clusters=5.3%, P=0.01, 519 Suppl. Table 6). 520 Black spruce populations showed a significant pattern of isolation by distance (IBD, Mantel’s 521 R=0.68, P=8×10-4, Suppl. Fig. 1). Genetic clusters were geographically separated in space (Fig. 522 1B), and the two main axes of spatial autocorrelation (Moran eigenvector maps, MEMs) 523 overlapped with the distribution of these clusters (Suppl. Fig. 2), agreeing with the IBD pattern, 524 and confirming that geographic distribution is an important contributor to genetic differentiation 525 among black spruce populations. The two main axes of spatial PCA, confirmed separation of 526 cluster West, and a split between clusters East and Central (Suppl. Fig. 3), lending support to 527 the presence of three genetic clusters. 528 Observed heterozygosity ranged between 0.078 and 0.101 per population (Fig. 1F). Most 529 populations from the West cluster showed lower heterozygosity and fewer polymorphic sites. 530 When heterozygosity was measured including only monomorphic sites within each population, 531 West cluster populations showed higher observed heterozygosity in spite of lower number of 532 polymorphic sites (Suppl. Fig. 4), showing that in West cluster heterozygosity is driven by a 533 large number of monomorphic sites rather than the higher fraction of homozygous sites. 534 Local adaptation of black spruce to climate 535 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 21 The three genetic clusters of black spruce occupied different climate niches. Population 536 locations were characterized in terms of climate by compiling 43 climate variables averaged 537 across 30 years (from 1961 to 1990, Suppl. Table 3, Suppl. Dataset 3). The first two principal 538 components (PCs) summarising those variables distinguished populations from the West, 539 Central, and East, with eastern locations being warmest and most humid and western ones 540 being driest and coldest (Fig. 2A). Partial RDA analysis conditioned on genetic structure, spatial 541 structure and climate variables showed that each class of variables explained only a small 542 portion of genetic variation on its own (1.9-6.7%, P<0.001, Suppl. Table 7). Climate explained 543 only 1.9% of variation on its own, with most of the variation overlapping with the spatial and 544 genetic structure, consistent with the patterns of isolation by environment. To identify which 545 environmental variables were best associated with distribution of genetic diversity, we run 546 Gradient Forest on population frequencies estimated for 1000 randomly selected variants. 547 Photoperiod, climate moisture index, precipitation, fall frost probability were among the most 548 important, besides the axes of spatial autocorrelation (MEMs) (Fig. 2A-B). 549 Measurements of multiple phenotypic traits in common gardens revealed that black spruce 550 populations from climates distant from the garden’s conditions tended to show reduced fitness 551 compared to those from similar climates, implying local adaptation to climate. Specifically, we 552 compared population means of 17 fitness traits measured in four common gardens (PR, ML, 553 CH, AC, see location on Fig. 1A, Suppl. Table 2), with the Euclidean climate distance of each 554 population from its origin (climate transfer distance). Height, biomass increment, and DBH 555 (diameter at breast height) showed the strongest negative correlations with climate distance, 556 consistent with local adaptation, although these correlations were consistently weaker in PR 557 (rho<-0.43, Fig. 2C-D). Survival rate and wood density were generally weakly correlated. Short-558 term biomass increments showed varying patterns depending on the age of the trees and the 559 common garden: more negative correlations were observed for later stages of growth in CH and 560 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 22 PR, and for younger stages in ML and AC (Fig. 2D). Also metrics of drought resilience varied, 561 with weakly negative correlations observed only fo r resilience (Rs) and resistance (Rt) in AC 562 and ML (Fig. 2D). 563 Testing genomic offset with common gardens 564 Because many phenotypic traits showed negative correlations with climate distance across 565 common gardens, we used the measured traits as ground truth to evaluate how accurately 566 genomic offset predicts fitness loss in black spruce, and determine the factors that influence 567 model performance (Fig. 3). We investigated the impact of subsequent data aspects on 568 genomic offset predictions, namely: 1) type and number of markers used for training the model, 569 2) choice of phenotypic trait and common garden for testing, 3) size of train and test datasets, 4) 570 species- vs cluster-based predictions, 5) choice of genetic clusters for training, and 6) choice of 571 climate variable sets for training. 572 Type and number of markers used for training the model 573 Because the genetic structure of black spruce is tightly linked with the distribution of climate, we 574 hypothesized that random markers selected from the genome and markers associated with 575 climate would generate similar genomic offsets. Indeed we found no systematic differences in 576 model performance when using subsets of random markers, all available markers, or markers 577 associated with three climate PCs, either corrected, or not-corrected for genetic structure (Fig. 578 4A, Suppl. Fig. 5). The only notable pattern was that using LFMM outliers consistently led to 579 more negative rho values in PR across nearly all traits. Genomic offsets generated with models 580 trained on different marker sets were very highly correlated, except for models trained on LFMM 581 markers, which were less correlated with other models across all gardens and traits (Suppl. 582 Figs. 6-7). LFMM outliers also showed lower heterozygosity than other types of markers (Suppl. 583 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 23 Fig. 8), suggesting that they are enriched for variants that are fixed in genetic clusters. Inclusion 584 of low-frequency markers did not have an effect on model performance (Fig. 4A). 585 Choice of phenotypic trait and common garden for testing 586 We found large differences in model performance among common gardens compared to 587 differences among traits (Fig. 4A, Suppl. Fig. 5). In ML, CH and AC, the stronger the trait was 588 correlated with climate transfer distance, the better it was correlated with genomic offset. The 589 best model performance was found in AC, in spite of the fact that the training dataset contained 590 the fewest populations (n=10), and most populations in the test dataset were not in the training 591 dataset. It was also the only garden where the model performed better in predicting height than 592 climate alone. On the contrary, in PR, genomic offsets were not correlated or they were 593 positively correlated with fitness (Fig. 4A, Suppl. Fig. 5). The average model performance on the 594 train data was often slightly better than that on the test data, suggesting that model performance 595 can decrease when projecting offsets for new trees. However the difference was very small and 596 for a few traits we saw the reverse, with test datasets leading to better model performances than 597 train datasets (Suppl. Fig. 5). This could reflect large variation of fitness values within 598 populations, and the fact that our “ground truth” depends on the sample size. 599 Size of train and test datasets 600 We tested how downsampling of either test or train datasets can impact model performance. 601 Given high fitness variation within populations, larger test or train samples (more trees per 602 population) should lead to more accurate estimates of population means, and subsequently 603 increase model performance. On the other hand, including more populations in the training 604 dataset should improve model performance by incorporating a wider range of climates into the 605 modeling of genotype-climate associations. We found that when trait population means in test 606 datasets were averaged across four or less individuals, average model performances usually 607 decreased in AC, ML and CH, but changes were very small, and mostly noticeable in AC 608 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 24 (Suppl. Fig. 9). Surprisingly average model performances increased in PR for smaller test 609 sample sizes (less positive rho, Suppl. Fig. 9). 610 Next, we investigated whether reducing the number of populations in the train dataset 611 decreased model performance. We found that in ML, CH and AC, downsampling to less than 10 612 populations usually decreased average model performance, with the strongest effect again in 613 AC (Suppl. Fig. 10). In PR, the effect was reversed, with less populations driving correlations 614 down (therefore improving performance), even though estimates of rho did not drop below zero. 615 Lastly, we tested the impact of the number of genotypes used to calculate population allele 616 frequencies (size of train population samples). Smaller train sample sizes had a minor negative 617 effect on average model performance in AC, they had a positive effect in PR, whereas the effect 618 in ML and CH was negligible or mixed (Suppl. Fig. 11). Overall, removing populations or 619 genotypes from train datasets tended to bring correlations of fitness-offset closer to zero. While 620 the impacts of removing genotypes were subtle, removing populations down to six could shift 621 the average model performance by 20%. Notably, adding more than ten populations to the 622 training set did not lead to any further improvements in average model performance across 623 most fitness traits. 624 Species- vs cluster-based predictions 625 We next investigated how the models predict fitness of individual genetic clusters. We found 626 large differences in model performance when clusters were tested separately (Fig. 4B, Suppl. 627 Fig. 12). Genomic offsets of the East cluster were negatively correlated with almost all traits in 628 each of three common gardens where they were tested. On the contrary, Central cluster was 629 correlated with fitness only in AC, and West cluster only in ML. Correlation between traits and 630 climate transfer distance showed similar patterns (Fig. 4B, Suppl. Figs. 13-14), although we still 631 observed discrepancy between climate distance and genomic offset correlations in PR. 632 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 25 Choice of genetic clusters for training 633 We investigated whether training Gradient Forest models independently on each cluster would 634 enhance the accuracy of cluster-specific fitness predictions and also enable cross-cluster 635 predictions. We found that the choice of a dataset for training (combined or separate clusters) 636 most of the time had very little effect on model performance with few exceptions (Fig. 5A). When 637 models were trained on a Central cluster, predictions of height in Central cluster slightly 638 improved in ML, and predictions of biomass increments in Central and West clusters slightly 639 improved in PR, although the changes were small and did not occur across all traits (Suppl. Fig. 640 15). Genomic offsets generated with models trained on different genetic clusters were positively 641 correlated but correlations varied across models, traits and gardens (Suppl. Fig. 16). For height 642 and biomass increment, the weakest correlations were observed in PR, between models trained 643 on the Central cluster and other models (Pearson correlation coefficient ranging between 0.61 644 and 0.86). 645 Choice of climate variable sets for training 646 A possible reason for inability to predict trait response of some of the clusters with genomic 647 offset, is that clusters may respond differently to different combinations of climate variables, and 648 these are not well-captured by the first three climate PCs. We therefore trained Gradient Forest 649 model with five variables (dPP, CMI, TP, WSTP, and fallMT) selected out of 15 (a subset of all 650 variables for which we had both past and future values and therefore could be used for future 651 offset predictions) that were identified as most important by the model. We found that in ML, 652 CH, and AC, the model performances dropped regardless of whether the models were trained 653 on combined datasets or single clusters only (Fig. 5B, Suppl. Fig. 17). On the contrary, in PR, all 654 models produced genomic offsets that were negatively correlated with the two traits in the West 655 cluster, and also with the biomass increment in the Central cluster, but not when fitness traits 656 were combined for the species. Most genomic offsets generated by models trained on clusters 657 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 26 or combined dataset were well correlated with each other (Pearson correlation coefficient >0.9), 658 except in CH, where correlation coefficients ranged between 0.63 and 0.99 (Suppl. Fig. 18). On 659 the contrary, when genomic offsets were compared between models trained on two separate 660 climate sets of variables, the correlation coefficients were much lower, ranging between -0.44 661 and 0.90 with a median as low as 0.50 (Suppl. Fig. 19). This suggests that different 662 combinations of climate variables may better capture climate black spruce cluster responses at 663 different locations. 664 Vulnerability of black spruce populations to future climate 665 Common garden experiments are ideal for testing hypotheses about the consequences of 666 assisted migration, but at the same time we are interested in predicting how natural populations 667 will perform in the face of changi ng conditions. To visualize the shift in climate across black 668 spruce distribution, we classified the climate state of each population in the past and across four 669 future periods (2022-40, 2041-60, 2061-80, and 2081-2100) under three carbon emission 670 scenarios into ten clusters (Suppl. Fig. 20). The West, Central, and East clusters were projected 671 to experience progressively warmer climates over time and along the south-north gradient. 672 (Suppl. Fig. 20). 673 We estimated black spruce genomic offset for the four time periods and a medium emission 674 scenario SSP2-4.5. Because our validation results showed discrepancy in performance 675 between western and central/eastern common gardens when using different climate variables, 676 we trained one model using West cluster populations and selected five climate variables to 677 estimate genomic offsets for the West cluster, and the second model using combined 678 populations and three climate PCs to estimate genomic offsets for the Central and East clusters. 679 In the western part of the black spruce distribution, genomic offsets were predicted to increase 680 across most of the West cluster range, reaching their highest values in the northwest (Fig. 6A). 681 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 27 In the central and eastern regions, the model indicated higher offsets in the central populations 682 relative to those in the east (Fig. 6B). 683

684 Distribution of genetic diversity in black spruce 685 Black spruce carries high levels of genetic variation and has large estimated population sizes 686 across its distribution due to its outcrossing mating system and wind pollination (Bouillé & 687 Bousquet, 2005; Isabel et al., 1995; Perry & Bousquet, 2001). Our results confirm that variation 688 is sorted among three widespread genetic clusters: West, Central, and East, with the West 689 cluster being most distinct (Gérardi et al., 2010; Jaramillo-Correa et al., 2004; Prunier et al., 690 2012). This structure is interpreted as the legacy of the last glacial maximum (Jaramillo-Correa 691 et al., 2004), although our data does not exclude the possibility that Central and East clusters 692 were derived from a single glacial refugium, and later diverged after migrating and adapting to 693 distinct environmental conditions. We also identified two isolated subclusters, one within Central 694 (WI) and one within East cluster (ME). While WIs range is most likely limited to a small area, in 695 case of ME it’s difficult to judge how widespread this subcluster is as it is situated at the 696 southern limits of the sampling range, but not at the limit of the species range. Pairwise FSTs 697 show these two subclusters are genetically differentiated from other clusters, from each other, 698 as well as from the red spruce, which suggests that they could be remnants of glacial 699 microrefugia that did not succeed in expanding across the landscape. Such cryptic refugia were 700 indeed identified in other North American species (Fernandez et al., 2021; Shafer et al., 2010). 701 We see some variation in the levels of observed heterozygosity with the East cluster 702 populations having higher, and the West cluster populations - lower heterozygosity. Even 703 though we applied measures to reduce the potential biases coming from low frequency variants, 704 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 28 missing data, and uneven sample sizes, heterozygosity estimates based on variant sites only 705 are driven in structured populations by the amount of divergence between those populations 706 which determines the amount of population-specific non-variant sites left in the dataset (Schmidt 707 et al., 2021). With this in mind, we also reported heterozygosity calculated only on the 708 polymorphic sites within each population, which shows that West cluster populations in spite of 709 fewer polymorphic sites, carry proportionally more heterozygous variants than other 710 populations. This would suggest that West populations do not necessarily have smaller effective 711 population sizes than East or Central populations. 712

to genomic offset in black spruce 713 Overall, our validation approach revealed that within common gardens, genomic offset 714 predictions were largely invariant to the choice of markers, the size of the training dataset, or the 715 genetic cluster on which model was trained. Genomic offset accuracies varied depending on the 716 common garden considered, fitness trait as well as the genetic cluster for which offset was 717 measured, and the selection of climate variables had the largest influence on both model 718 performance and the ranking of population offsets. 719 Variation in genomic offset performance across gardens and clusters 720 One of the main assumptions of genomic offset metrics is that local populations are under 721 stabilizing selection for optimal phenotypes (Gain et al., 2023). With common garden 722 experiments we demonstrated that populations are most likely adapted to their local climate 723 conditions, however population response to the same set of climate variables varied depending 724 on the garden, phenotypic trait, and genetic cluster to which population belonged. In other 725 words, two populations can respond differently to increased humidity, depending on their 726 genetic background, and whether they grow in cold sites with a short growing season or warmer 727 sites with long-growing seasons. This is also reflected in estimated genomic offsets, which most 728 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 29 of the time are correlated with fitness traits to the same extent as climate transfer distance. 729 Notably, in contrast to previous validation studies (Fitzpatrick et al., 2021; Lachmuth et al., 730 2024) we rarely observed genomic offset outperforming climate transfer distance, with an 731 exception of the AC garden, where fitness-offset correlations reached -0.8 and were better than 732 correlation coefficients based on climate transfer distance alone. Genomic offset is expected to 733 model fitness loss better than climate, because climate variables are weighted by their effects 734 on association with genetic variation (Gain et al., 2023). However when selected climate 735 variables are known to be associated with adaptive variation, climate distances can actually 736 better predict fitness than genomic offsets (Láruson et al., 2022). We calculated climate transfer 737 distance for a composite of >40 climate variables, specifically because fitness traits were better 738 correlated with distance calculated in this way rather than using single climate variables. 739 Consequently, subsequent Gradient Forest models were run using three major principal 740 components (climate PCs) derived from a large set of climate variables, instead of a set of 741 separate climate variables. The rationale was that even though climate PCs are harder to 742 interpret, they better capture local climate niches of black spruce. Such choice of variables 743 could explain why we obtained usually similar results using climate and genomic offset. 744 Similarly, in Douglas-fir, models trained on climate variables identified as influential for fitness in 745 common garden experiments proved more accurate than those trained on other variables (Lind 746 et al., 2024). 747 In some cases, genomic offset showed poorer performance than climate, most notably for PR, 748 in which genomic offsets were either not correlated or positively correlated with fitness traits. 749 This positive association in PR seems to partially result from a statistical artifact, where the 750 correlation is driven by unexpected mean fitness differences between West and Central 751 clusters. Even though both clusters show negative correlation with climate distance, trees from 752 the Central cluster grow on average better than those from the West cluster (Suppl. Figs. 13-14, 753 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 30 (Robert et al., 2024)), and genomic offsets are also higher for the Central cluster, resulting in 754 species-level correlations being positive (Suppl. Fig. 21). Presence of these noncrossover 755 genotype-by-environment interactions could also explain the positive fitness-offset correlations 756 we occasionally observed within clusters. For instance, in PR, if hybrid individuals display 757 intermediate traits between West and Central clusters, their inclusion in the West cluster might 758 have artificially inflated both fitness estimates and genomic offsets. Additional populations from 759 the western part of Canada would be helpful in testing this hypothesis. Another, non-mutually 760 exclusive explanation is that in western Canada, other sets of climate factors play a role in local 761 adaptation (see section “ The choice of climate variables for training impacts model 762 performance”). Soil conditions and microbiota can also impact host phenotype and how species 763 respond to environmental stressors and adapt to new environments (Allsup, George, and 764 Lankau 2023; Henry et al. 2021), therefore including these factors in the future can potentially 765 improve trait response curves. Overall, these results suggest that even when genomic offset 766 can accurately estimate maladaptation of populations, such maladaptation can be manifested 767 differently in terms of phenotype across the species landscape, especially in case of species 768 with very wide distributions, strong genetic structure and potential morphologically distinct 769 ecotypes, therefore it will not predict correctly a specific fitness-related phenotype. This result 770 resonates with recent findings conducted in genetically structured maritime pine ( Pinus pinaster 771 Ait.), which found mixed results of fitness-offset correlations across multiple common gardens 772 (Archambeau et al., 2025). 773 Variation in genomic offset performance across fitness traits 774 Model performances varied across traits, which in our case simply reflects the variation in trait 775 responses to climate. Height, total biomass and DBH in general showed the strongest negative 776 correlations with climate and offsets. These functional traits are often used as a proxy for fitness 777 in tree species, because they respond to climate in common garden experiments and show 778 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 31 moderate heritability (Fitzpatrick et al. 2021; Housset et al. 2018; Girardin et al. 2021; Robert et 779 al. 2024). Interestingly, short-term biomass increments showed contrasting performances 780 across the gardens depending on the tree age. In the two warmest gardens (ML and AC), early 781 intervals of biomass increment were better associated with climate distance (and genomic 782 offsets), whereas in the two coldest ones (PR and CH), the later intervals showed the strongest 783 associations, both at the species and cluster level. Low correlations between juvenile and 784 mature traits have been previously observed in conifers (Depardieu et al., 2021). In terms of 785 genomic offset, these results only partially agree with recent findings in red spruce, which 786 showed stronger associations for juvenile-stage phenotypes, when both juvenile- and late-stage 787 phenotypes were compared (Verrico et al., 2025). As red spruce gardens were located at lower 788 latitudes than those of black spruce, it could suggest that temperature or photoperiod have an 789 impact on early vs late-stage response of trees. We did not use survival in genomic offset 790 validations, since we found it responded weakly to climate, with the strongest response in PR 791 and AC. While survival can be a useful metric of fitness (Archambeau et al., 2025; Lind et al., 792 2024), in this case it was measured only for seeds that managed to germinate and develop into 793 seedlings, therefore it did not incorporate early-stage mortality which could have influenced the 794 overall survival estimates. 795 While mean drought resilience metrics varied across gardens and clusters, they were not 796 strongly associated with climate distance or genomic offsets, with the strongest negative 797 associations observed for Rt (resistance) and Rs (resilience) in ML. This appears to be 798 consistent with independent models of recovery to resistance trade-offs (the trade-off between 799 resistance to drought event and ability to recover afterwards, Robert et al. in prep.), which show 800 that in general the effects of population climate or genetic ancestries were weak. Nevertheless, 801 in ML, the population variation in trade-off was impacted by the growth degree days. 802 Specifically, populations coming from regions with shorter growing seasons (more distant from 803 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 32 the garden) were associated with higher-resistance-lower-recovery strategy. Given that most 804 populations planted in ML come from regions with shorter growing periods, this translates to 805 negative association between resistance and climate distance. 806 The number or type of markers has little impact on model performance 807 We did not find any consistent patterns in the choice of markers on offset predictions, neither in 808 terms of the marker number, distribution of marker allele frequencies or the strength of 809 association with climate variables. This is generally consistent with many studies validating 810 genomic offset using Gradient Forest in common gardens or via simulations (Fitzpatrick et al., 811 2021; Láruson et al., 2022; Lind et al., 2024; Lind & Lotterhos, 2025a; Verrico et al., 2025). The 812 likely explanation is that much of segregating variation in black spruce is linked to the 813 environment as a consequence of spatial distribution of genetic clusters, and patterns of 814 isolation by distance driven by expansion of these clusters from glacial refugia through the 815 landscape. Such patterns of isolation by environment would result in most genomic variants 816 being linked to climate via monotonic associations, with climate-adapted loci providing no 817 additional advantage. Even though, intuitively larger number of makers should lead to 818 identification of more and stronger associations with climate, this does not result in 819 improvements of genomic offsets, as has been demonstrated using simulations (Lind & 820 Lotterhos, 2025a). 821 Most variation in genomic offset performance was observed when using LFMM outliers, and in 822 particular in PR where fitness-offset correlations tended to decrease for these markers in almost 823 all traits. LFMM outliers showed lower heterozygosity compared to other markers, and after 824 inspection of outlier frequencies, most of them captured fixed or nearly fixed differences 825 between East and Central/West clusters. As a consequence, in PR, which comprises only West 826 and Central cluster populations, LFMM outliers did not differentiate well between the two 827 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 33 clusters, leading to overlapping ranges of genomic offsets for both clusters without changing 828 their fitness values (Suppl. Fig. 21). 829 The size of training and testing datasets has little impact on model performance 830 We found that size limitations on training and test populations showed a very minor effect on 831 genomic offset predictions. Across three common gardens (CH, ML, AC), typically 832 downsampling to less than four to six genotypes per population led to slightly lower model 833 performances. Performances also usually dropped when sizes of testing populations were very 834 small (two or four samples). Increased performance in genomic offset with increasing test or 835 train population sizes is expected in cases when populations are genetically and phenotypically 836 diverse, because averaging across larger populations reduces the error caused by 837 measurement or other non-genetic effects. Indeed in black spruce, a large proportion of genetic 838 and phenotypic variance is distributed within populations (Rossi & Bousquet, 2014). This 839 observation reflects simulation results showing that genomic offsets based on population fitness 840 means tend to have higher performance compared to those based on individual fitness 841 measurements (Lind & Lotterhos, 2025a). 842 Increasing the number of populations in the training dataset can improve the modeling of 843 genotype-environment associations (Láruson et al., 2022). Surprisingly, we found a very small 844 impact of reducing the number of populations in the training dataset below ten, with the 845 strongest impact being in AC. Similarly, training models on populations coming only from a 846 single genetic cluster did not largely affect how genomic offsets were correlated with fitness 847 across species but also across populations from other clusters. These different models also 848 produced genomic offsets that were strongly correlated with each other, indicating similar spatial 849 distribution patterns and ranking. One interpretation of this is that given high genetic diversity of 850 black spruce in general and within populations, even a small number of randomly sampled 851 populations capture enough genetic variation and their associations with climate necessary to 852 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 34 train the model. It could suggest that linear approximations of Gradient Forest to novel locations 853 and climates prove sufficient in the case of black spruce. 854 Interestingly, all our sensitivity tests revealed an opposite pattern in PR. Specifically, reducing 855 the population sizes of the training and test datasets, as well as the number of populations in 856 the training set, generally decreased the otherwise positive fitness-offset correlations observed 857 across most traits. As discussed above, such positive correlations likely arise from cluster-858 specific differences in mean fitness, which may be expressed in hybrid individuals. If this 859 explanation holds, reducing sample sizes could remove hybrids and thereby eliminate these 860 spurious correlations. 861 The choice of climate variables for training impacts model performance 862 The choice of climate variables used to train the models had the largest influence on fitness-863 offset correlations, and on how genomic offsets were ranked spatially. Models trained on three 864 climate PCs were best at predicting fitness of the East cluster and the whole species fitness in 865 eastern common gardens (CH, ML and AC), whereas models using selected five climate 866 variables were better at predicting fitness of the West cluster in western common garden (PR). 867 Temperature and precipitation interact nonlinearly to influence black spruce growth 868 (D’Orangeville et al., 2018). In the east, where precipitation is hi gh, low temperatures in the 869 north are the major limitation of black spruce growth, while in the south, high temperatures can 870 limit growth by higher evapotr anspiration (D’Orangeville et al ., 2018; Girardin et al., 2016; 871 Lesven et al., 2024). In the western and central parts of Canada, precipitation is much lower, 872 and moisture availability is the major determinant of growth (Lesven et al., 2024; Price et al., 873 2013). This appears to pose a challenge to genomic offset models which do not take into 874 account interactions between climate variables across the landscape. As a consequence we 875 were not able to provide one model to accurately predict population responses to climate across 876 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 35 the four common gardens. This also underlines the importance of finding the most influential 877 environmental variables for genomic offset estimation. 878 Future predictions in black spruce 879 Common garden experiments are ideal for testing tree responses to assisted migration (“local-880 foreign” fitness offset sensu Lotterhos, 2024). But in the case of black spruce we are also 881 interested in predicting vulnerability of natural populations across the landscape to future 882 climate changes (“home-away” fitness offset sensu Lotterhos, 2024). In the context of common 883 garden experiments, validation of the second scenario would require fitness-offset comparisons 884 between populations planted in their native climate and the same populations planted in distant 885 climates, preferentially in climates close to those anticipated in the future. In absence of these, 886 one can attempt to leverage information derived from garden-specific models developed here to 887 estimate genomic offsets to future climate. Three key considerations should be kept in mind 888 when applying this approach. First, this and previous studies (Robert et al., 2024) identified 889 asymmetry in reciprocal transplants of populations between western and eastern parts of black 890 spruce distribution. Populations from the West cluster consistently exhibit lower growth across 891 all gardens, including those closest to their origin, violating the validation requirement that local 892 populations should exhibit the highest relative fitness. This renders genomic offsets 893 uninformative in terms of predicting fitness traits at the species level. Second, as mentioned 894 above, populations from the western part of the species range (represented by the PR common 895 garden) appear to respond to a different set of climate variables as those from central or eastern 896 Canada (common gardens CH, ML and AC). This argues for using separate models focusing 897 either on the Central/East or on the West cluster of black spruce, with each model using a 898 different set of environmental variables for training. Finally, we found some variation in climate 899 response across individual clusters measured within common gardens. Even between the two 900 most related clusters, Central and East, the Central cluster exhibited consistently weaker 901 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 36 fitness-climate and fitness-offset correlations on average. More investigations need to be done 902 to determine whether such differences are driven by sampling differences, hybridization 903 between clusters, between species (hybridization between the East cluster and red spruce), or 904 different responses of the clusters to climate. 905 The two projections of genomic offsets to future climate in black spruce largely agreed with the 906 literature. The first model (model Central/East) identified higher genomic offsets in the Central 907 compared to the East cluster. This is consistent with previous studies indicating that western 908 and central black spruce populations are expected to experience the greatest mortality due to 909 climate change (Lesven et al., 2024), whereas eastern populations are projected to be less 910 affected and to respond with a delay (D’Orangeville et al., 2018; Girardin et al., 2021; Robert et 911 al., 2024). The second model (model West) suggested a nonuniform distribution of genomic 912 offsets, ultimately leading to the highest offsets for north-west populations, which is supported 913 by expectations of population declines in that region (Walker et al., 2015). Importantly, genomic 914 offset methods can predict only losses in fitness, not potential gains. Previous models have 915 suggested that productivity at northernmost black spruce locations may temporarily increase 916 with rising temperatures (D’Orangeville et al., 2018; Robert et al., 2024). However, such 917 changes are expected to be short-lived. Therefore, genomic offset estimates should be 918 interpreted with caution. In summary, although we are unable to compare genomic offsets 919 across the entire black spruce distribution, the Gradient Forest models provided projections that 920 at least, on a broad scale, aligned with independent evidence. 921 922

923 By sequencing and phenotyping distribution-wide populations of black spruce in four long-term 924 common garden experiments we validated genomic offsets as predictors of fitness in this 925 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 37 genetically structured species. Gradient Forest models appeared to be constrained in predicting 926 climate responses across the entire distribution of black spruce. This limitation arises, first, 927 because similar ranks in cluster relative fitness are retained between gardens, which prevents 928 reliable predictions of fitness across the species, and second, because species responses to 929 climate vary across the landscape not only as a function of individual climatic variables but also 930 of their interactions. On the other hand, models trained on a specific set of climate variables 931 proved to be very robust to the choice of markers, but also to the amount and range of 932 populations, as well as number of genotypes per population used in model training. Despite the 933 uncertainty associated with testing only four independent sites (and thus only four 934 environments), we identified two models suitable for predicting genomic offsets in the western 935 and central/eastern parts of the black spruce distribution, respectively. Future directions should 936 include establishment of new common gardens, especially in the west, where population 937 responses to climate are less understood. Moreover, good quality individual genotypes, and 938 denser sampling could prove useful in improving genomic offset predictions on a more local 939 scale, through more accurate modeling of fine-scale genotype-environment associations. 940 Data Archiving Statement 941 Data for this study are available at: to be completed after manuscript is accepted for publication. 942 943 Literature cited 944 Aguirre-Liguori, J. A., Ramírez-Barahona, S., & Gaut, B. S. (2021). The evolutionary genomics 945 of species’ responses to climate change. Nature Ecology & Evolution, 5(10), 1350–1360. 946 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 38 Alexander, D. H., Novembre, J., & Lange, K. (2009). Fast model-based estimation of ancestry in 947 unrelated individuals. Genome Research, 19(9), 1655–1664. 948 Archambeau, J., Benito Garzón, M., de-Miguel, M., Changenet, A., Bagnoli, F., Barraquand, F., 949 Marchi, M., Vendramin, G. G., Cavers, S., Perry, A., & González-Martínez, S. C. (2025). 950 Evaluating genomic offset predictions in a forest tree with high population genetic structure. 951 The American Naturalist, 739045. https://doi.org/10.1086/739045 952 Aubin, I., Boisvert-Marsh, L., Kebli, H., McKenney, D., Pedlar, J., Lawrence, K., Hogg, E. H., 953 Boulanger, Y., Gauthier, S., & Ste-Marie, C. (2018). Tree vulnerability to climate change: 954 improving exposure-based assessments using traits as indicators of sensitivity. Ecosphere 955 (Washington, D.C), 9(2), e02108. 956 Baltzer, J. L., Day, N. J., Walker, X. J., Greene, D., Mack, M. C., Alexander, H. D., Arseneault, 957 D., Barnes, J., Bergeron, Y., Boucher, Y., Bourgeau-Chavez, L., Brown, C. D., Carrière, S., 958 Howard, B. K., Gauthier, S., Parisien, M.-A., Reid, K. A., Rogers, B. M., Roland, C., … 959 Johnstone, J. F. (2021). Increasing fire and the decline of fire adapted black spruce in the 960 boreal forest. Proceedings of the National Academy of Sciences of the United States of 961 America, 118(45), e2024872118. 962 Bellemin-Noël, B., Bourassa, S., Despland, E., De Grandpré, L., & Pureswaran, D. S. (2021). 963 Improved performance of the eastern spruce budworm on black spruce as warming 964 temperatures disrupt phenological defences. Global Change Biology, 27(14), 3358–3366. 965 Bergmann, F. (1978). The allelic distribution at an acid phosphatase locus in Norway spruce 966 (Picea abies) along similar climatic gradients. Theoretical and Applied Genetics, 52(2), 57–967 64. 968 Bivand, R. S., & Wong, D. W. S. (2018). Comparing implementations of global and local 969 indicators of spatial association. Test (Madrid, Spain), 27(3), 716–748. 970 Bouillé, M., & Bousquet, J. (2005). Trans-species shared polymorphisms at orthologous nuclear 971 gene loci among distant species in the conifer Picea (Pinaceae): implications for the long-972 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 39 term maintenance of genetic diversity in trees. American Journal of Botany, 92(1), 63–73. 973 Bradley St Clair, J., & Howe, G. T. (2007). Genetic maladaptation of coastal Douglas/i1fir 974 seedlings to future climates. Global Change Biology, 13(7), 1441–1454. 975 Brandt, J. P., Flannigan, M. D., Maynard, D. G., Thompson, I. D., & Volney, W. J. A. (2013). An 976

to Canada’s boreal zone: ecosystem processes, health, sustainability, and 977 environmental issues. Environmental Reviews, 21(4), 207–226. 978 Burns, R. M., & Honkala, B. H. (1990). Silvics of North America: Conifers (Forest Service, 979 United States Department of Agriculture (ed.)). Oxford. 980 Capblancq, T., & Forester, B. R. (2021). Redundancy analysis: A Swiss Army Knife for 981 landscape genomics. Methods in Ecology and Evolution, 12(12), 2298–2309. 982 Capblancq, T., Luu, K., Blum, M. G. B., & Bazin, E. (2018). Evaluation of redundancy analysis to 983 identify signatures of local adaptation. Molecular Ecology Resources, 18(6), 1223–1233. 984 Capblancq, T., Tauzin, A., Vigouroux, Y., Cubry, P., & François, O. (2025). Linking genomic 985 offset statistics to the shape of selection gradients. The American Naturalist, 739079. 986 https://doi.org/10.1086/739079 987 Depardieu, C., Gérardi, S., Nadeau, S., Parent, G. J., Mackay, J., Lenz, P., Lamothe, M., 988 Girardin, M. P., Bousquet, J., & Isabel, N. (2021). Connecting tree-ring phenotypes, genetic 989 associations and transcriptomics to decipher the genomic architecture of drought 990 adaptation in a widespread conifer. Molecular Ecology, 30(16), 3898–3917. 991 Depardieu, C., Girardin, M. P., Nadeau, S., Lenz, P., Bousquet, J., & Isabel, N. (2020). Adaptive 992 genetic variation to drought in a widely distributed conifer suggests a potential for 993 increasing forest resilience in a drying climate. The New Phytologist, 227(2), 427–439. 994 D’Orangeville, L., Houle, D., Duchesne, L., Phillips, R. P., Bergeron, Y., & Kneeshaw, D. (2018). 995 Beneficial effects of climate warming on boreal tree growth may be transitory. Nature 996 Communications, 9(1), 3213. 997 Ellis, N., Smith, S. J., & Pitcher, C. R. (2012). Gradient forests: calculating importance gradients 998 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 40 on physical predictors. Ecology, 93(1), 156–168. 999 Evans, L. M., Slavov, G. T., Rodgers-Melnick, E., Martin, J., Ranjan, P., Muchero, W., Brunner, 1000 A. M., Schackwitz, W., Gunter, L., Chen, J.-G., Tuskan, G. A., & DiFazio, S. P. (2014). 1001 Population genomics of Populus trichocarpa identifies signatures of selection and adaptive 1002 trait associations. Nature Genetics, 46(10), 1089–1096. 1003 Fernandez, M. C., Hu, F. S., Gavin, D. G., Lafontaine, G., Heath, K. D., & Vanderpoorten, A. 1004 (2021). A tale of two conifers: Migration across a dispersal barrier outpaced regional 1005 expansion from refugia. Journal of Biogeography, 48(9), 2133–2143. 1006 Fitzpatrick, M. C., Chhatre, V. E., Soolanayakanahally, R. Y., & Keller, S. R. (2021). 1007 Experimental support for genomic prediction of climate maladaptation using the machine 1008 learning approach Gradient Forests. Molecular Ecology Resources, 21(8), 2749–2765. 1009 Fitzpatrick, M. C., & Keller, S. R. (2015). Ecological genomics meets community-level modelling 1010 of biodiversity: mapping the genomic landscape of current and future environmental 1011 adaptation. Ecology Letters, 18(1), 1–16. 1012 Fitzpatrick, M. C., Keller, S. R., & Lotterhos, K. E. (2025). The challenge of genomic forecasting 1013 in an era of global change. The American Naturalist, 738891. 1014 https://doi.org/10.1086/738891 1015 Fortin, M., Lavoie, J.-F., Régnière, J., & Saint-Amant, R. (2022). A Web API for weather 1016 generation and pest development simulation in North America. Environmental Modelling & 1017 Software: With Environment Data News, 157(105476), 105476. 1018 Gain, C., & François, O. (2021). LEA 3: Factor models in population genetics and ecological 1019 genomics with R. Molecular Ecology Resources, 21(8), 2738–2748. 1020 Gain, C., Rhoné, B., Cubry, P., Salazar, I., Forbes, F., Vigouroux, Y., Jay, F., & François, O. 1021 (2023). A quantitative theory for genomic offset statistics. Molecular Biology and Evolution, 1022 40(6), msad140. 1023 Gauthier, S., Bernier, P., Kuuluvainen, T., Shvidenko, A. Z., & Schepaschenko, D. G. (2015). 1024 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 41 Boreal forest health and global change. Science (New York, N.Y.), 349(6250), 819–822. 1025 Gérardi, S., Jaramillo-Correa, J. P., Beaulieu, J., & Bousquet, J. (2010). From glacial refugia to 1026 modern populations: new assemblages of organelle genomes generated by differential 1027 cytoplasmic gene flow in transcontinental black spruce: ASSEMBLAGES OF ORGANELLE 1028 GENOMES. Molecular Ecology, 19(23), 5265–5280. 1029 Girardin, M. P., Guo, X. J., Marchand, W., & Depardieu, C. (2024). Unravelling the 1030 biogeographic determinants of tree growth sensitivity to freeze and drought in Canada’s 1031 forests. The Journal of Ecology, 112(4), 848–869. 1032 Girardin, M. P., Hogg, E. H., Bernier, P. Y., Kurz, W. A., Guo, X. J., & Cyr, G. (2016). Negative 1033 impacts of high temperatures on growth of black spruce forests intensify with the 1034 anticipated climate warming. Global Change Biology, 22(2), 627–643. 1035 Girardin, M. P., Isabel, N., Guo, X. J., Lamothe, M., Duchesne, I., & Lenz, P. (2021). Annual 1036 aboveground carbon uptake enhancements from assisted gene flow in boreal black spruce 1037 forests are not long-lasting. Nature Communications, 12(1), 1169. 1038 Global Forest Resources Assessment 2020. (n.d.). GlobalFRA. Retrieved September 21, 2025, 1039 from https://www.fao.org/forest-resources-assessment/past-assessments/fra-2020/en 1040 Gosselin, T. (2020). thierrygosselin/radiator: update. Zenodo. 1041 https://doi.org/10.5281/ZENODO.3687060 1042 Gosselin, T., Anderson, E. C., Ferchaud, A.-L., & Bar, I. (2019). thierrygosselin/assigner: 1043 assigner 0.5.6 2019-05-01. Zenodo. https://doi.org/10.5281/ZENODO.2656728 1044 Guénard, G., & Legendre, P. (2022). Hierarchical clustering with contiguity constraint in R. 1045 Journal of Statistical Software, 103(7), 1–26. 1046 Guo, X., Klisz, M., Puchałka, R., Silvestro, R., Faubert, P., Belien, E., Huang, J., & Rossi, S. 1047 (2022). Common/i1garden experiment reveals clinal trends of bud phenology in black 1048 spruce populations from a latitudinal gradient in the boreal forest. The Journal of Ecology, 1049 110(5), 1043–1053. 1050 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 42 Hogg, E. H., Barr, A. G., & Black, T. A. (2013). A simple soil moisture index for representing 1051 multi-year drought impacts on aspen productivity in the western Canadian interior. 1052 Agricultural and Forest Meteorology, 178-179, 173–182. 1053 Housset, J. M., Nadeau, S., Isabel, N., Depardieu, C., Duchesne, I., Lenz, P., & Girardin, M. P. 1054 (2018). Tree rings provide a new class of phenotypes for genetic associations that foster 1055 insights into adaptation of conifers to climate change. The New Phytologist, 218(2), 630–1056 645. 1057 IPCC. (2023). Climate Change 2023: Synthesis Report. Contribution of Working Groups I, II and 1058 III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change 1059 (Core Writing Team, H. Lee and J. Romero (Ed.); pp. 35–115). IPCC. 1060 Isaac-Renton, M. G., Roberts, D. R., Hamann, A., & Spiecker, H. (2014). Douglas-fir plantations 1061 in Europe: a retrospective test of assisted migration to address climate change. Global 1062 Change Biology, 20(8), 2607–2617. 1063 Isabel, N., Beaulieu, J., & Bousquet, J. (1995). Complete congruence between gene diversity 1064 estimates derived from genotypic data at enzyme and random amplified polymorphic DNA 1065 loci in black spruce. Proceedings of the National Academy of Sciences of the United States 1066 of America, 92(14), 6369–6373. 1067 Jaramillo-Correa, J. P., Beaulieu, J., & Bousquet, J. (2004). Variation in mitochondrial DNA 1068 reveals multiple distant glacial refugia in black spruce (Picea mariana), a transcontinental 1069 North American conifer. Molecular Ecology, 13(9), 2735–2747. 1070 Jombart, T. (2008). adegenet: a R package for the multivariate analysis of genetic markers. 1071 Bioinformatics , 24(11), 1403–1405. 1072 Kamvar, Z. N., Tabima, J. F., & Grünwald, N. J. (2014). Poppr: an R package for genetic 1073 analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ, 2, 1074 e281. 1075 Kilian, A., Wenzl, P., Huttner, E., Carling, J., Xia, L., Blois, H., Caig, V., Heller-Uszynska, K., 1076 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 43 Jaccoud, D., Hopper, C., Aschenbrenner-Kilian, M., Evers, M., Peng, K., Cayla, C., Hok, P., 1077 & Uszynski, G. (2012). Diversity arrays technology: a generic genome profiling technology 1078 on open platforms. Methods in Molecular Biology (Clifton, N.J.), 888, 67–89. 1079 Lachmuth, S., Capblancq, T., Prakash, A., Keller, S. R., & Fitzpatrick, M. C. (2024). Novel 1080 genomic offset metrics integrate local adaptation into habitat suitability forecasts and inform 1081 assisted migration. Ecological Monographs, 94(1), e1593. 1082 Láruson, Á. J., Fitzpatrick, M. C., Keller, S. R., Haller, B. C., & Lotterhos, K. E. (2022). Seeing 1083 the forest for the trees: Assessing genetic offset predictions from gradient forest. 1084 Evolutionary Applications, 15(3), 403–416. 1085 Lavoie, J., Bourgault, P., Smith, T. J., Logan, T., Leduc, M., Caron, L.-P., Gammon, S., & Braun, 1086 M. (2024). An ensemble of bias-adjusted CMIP6 climate simulations based on a high-1087 resolution North American reanalysis. Scientific Data, 11(1), 64. 1088 Lesven, J. A., Druguet Dayras, M., Cazabonne, J., Gillet, F., Arsenault, A., Rius, D., & 1089 Bergeron, Y. (2024). Future impacts of climate change on black spruce growth and 1090 mortality: review and challenges. Environmental Reviews, 32(2), 214–230. 1091 Lind, B. M., Candido-Ribeiro, R., Singh, P., Lu, M., Obreht Vidakovic, D., Booker, T. R., 1092 Whitlock, M. C., Yeaman, S., Isabel, N., & Aitken, S. N. (2024). How useful are genomic 1093 data for predicting maladaptation to future climate? Global Change Biology, 30(4), e17227. 1094 Lind, B. M., & Lotterhos, K. E. (2025a). A comparison of genomic forecasts based on genotypes 1095 versus allele frequencies. The American Naturalist, 739098. https://doi.org/10.1086/739098 1096 Lind, B. M., & Lotterhos, K. E. (2025b). The accuracy of predicting maladaptation to new 1097 environments with genomic data. Molecular Ecology Resources, 25(4), e14008. 1098 Liu, Q., Peng, C., Schneider, R., Cyr, D., McDowell, N. G., & Kneeshaw, D. (2023). Drought-1099 induced increase in tree mortality and corresponding decrease in the carbon sink capacity 1100 of Canada’s boreal forests from 1970 to 2020. Global Change Biology, 29(8), 2274–2285. 1101 Lloret, F., Keeling, E. G., & Sala, A. (2011). Components of tree resilience: effects of successive 1102 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 44 low-growth episodes in old ponderosa pine forests. Oikos (Copenhagen, Denmark), 1103 120(12), 1909–1920. 1104 Lortie, C. J., & Hierro, J. L. (2022). A synthesis of local adaptation to climate through reciprocal 1105 common gardens. The Journal of Ecology, 110(5), 1015–1021. 1106 Lotterhos, K. E. (2024). Interpretation issues with “genomic vulnerability” arise from conceptual 1107 issues in local adaptation and maladaptation. Evolution Letters, 8(3), 331–339. 1108 Marchand, W., Depardieu, C., Campbell, E. M., Bousquet, J., & Girardin, M. P. (2025). Long-1109 term temporal divergence in post-drought resilience decline between deciduous and 1110 evergreen tree species. Global Change Biology, 31(7), e70330. 1111 Marchand, W., Girardin, M. P., Hartmann, H., Lévesque, M., Gauthier, S., & Bergeron, Y. 1112 (2021). Contrasting life-history traits of black spruce and jack pine influence their 1113 physiological response to drought and growth recovery in northeastern boreal Canada. The 1114 Science of the Total Environment, 794(148514), 148514. 1115 Marquis, B., Bergeron, Y., Houle, D., Leduc, M., & Rossi, S. (2022). Variability in frost 1116 occurrence under climate change and consequent risk of damage to trees of western 1117 Quebec, Canada. Scientific Reports, 12(1), 7220. 1118 Mijangos, J. L., Gruber, B., Berry, O., Pacioni, C., & Georges, A. (2022). dartR v2: An 1119 accessible genetic analysis platform for conservation, ecology and agriculture. Methods in 1120 Ecology and Evolution / British Ecological Society. https://doi.org/10.1111/2041-210x.13918 1121 Morgenstern, E. K. (1978). Range-wide genetic variation of black spruce. Canadian Journal of 1122 Forest Research, 8(4), 463–473. 1123 Mura, C., Buttò, V., Silvestro, R., Deslauriers, A., Charrier, G., Raymond, P., & Rossi, S. (2022). 1124 The early bud gets the cold: Diverging spring phenology drives exposure to late frost in a 1125 Picea mariana [(Mill.) BSP] common garden. Physiologia Plantarum, 174(6), e13798. 1126 Oksanen, J., Simpson, G., Blanchet, F., Kindt, R., Legendre, P., Minchin, P., O’hara, R., 1127 Solymos, P., Stevens, M., Szoecs, E., Wagner, H., Barbour, M., Bedward, M., Bolker, B., 1128 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 45 Borcard, D., Borman, T., Carvalho, G., Chirico, M., De Caceres, M., … Weedon, J. (2025). 1129 vegan: Community Ecology Package (Version R package version 2.8-0). 1130 https://vegandevs.github.io/vegan/ 1131 Oleksyn, J., Modrzýnski, J., Tjoelker, M. G., Zytkowiak, R., Reich, P. B., & Karolewski, P. 1132 (1998). Growth and physiology of Picea abies populations from elevational transects: 1133 common garden evidence for altitudinal ecotypes and cold adaptation. Functional Ecology, 1134 12(4), 573–590. 1135 Patsiou, T. S., Shestakova, T. A., Klein, T., di Matteo, G., Sbay, H., Chambel, M. R., Zas, R., & 1136 Voltas, J. (2020). Intraspecific responses to climate reveal nonintuitive warming impacts on 1137 a widespread thermophilic conifer. The New Phytologist, 228(2), 525–540. 1138 Pebesma, E. (2018). Simple features for R: Standardized support for spatial vector data. The R 1139 Journal, 10(1), 439. 1140 Perry, D. J., & Bousquet, J. (2001). Genetic diversity and mating system of post-fire and post-1141 harvest black spruce: an investigation using codominant sequence-tagged-site (STS) 1142 markers. Canadian Journal of Forest Research, 31(1), 32–40. 1143 Prasad, A., Pedlar, J., Peters, M., McKenney, D., Iverson, L., Matthews, S., & Adams, B. (2020). 1144 Combining US and Canadian forest inventories to assess habitat suitability and migration 1145 potential of 25 tree species under climate change. Diversity & Distributions, 26(9), 1142–1146 1159. 1147 Price, D. T., Alfaro, R. I., Brown, K. J., Flannigan, M. D., Fleming, R. A., Hogg, E. H., Girardin, 1148 M. P., Lakusta, T., Johnston, M., McKenney, D. W., Pedlar, J. H., Stratton, T., Sturrock, R. 1149 N., Thompson, I. D., Trofymow, J. A., & Venier, L. A. (2013). Anticipating the consequences 1150 of climate change for Canada’s boreal forest ecosystems. Environmental Reviews, 21(4), 1151 322–365. 1152 Prunier, J., Gérardi, S., Laroche, J., Beaulieu, J., & Bousquet, J. (2012). Parallel and lineage-1153 specific molecular adaptation to climate in boreal black spruce: LINEAGE-SPECIFIC 1154 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 46 ADAPTIVE SNP IN BLACK SPRUCE. Molecular Ecology, 21(17), 4270–4286. 1155 R Core Team. (2021). R: A language and environment for statistical computing. R Foundation 1156 for Statistical Computing. https://www.R-project.org/ 1157 Rellstab, C., Dauphin, B., & Exposito-Alonso, M. (2021). Prospects and limitations of genomic 1158 offset in conservation management. Evolutionary Applications, 14(5), 1202–1212. 1159 Rellstab, C., Zoller, S., Walthert, L., Lesur, I., Pluess, A. R., Graf, R., Bodénès, C., Sperisen, C., 1160 Kremer, A., & Gugerli, F. (2016). Signatures of local adaptation in candidate genes of oaks 1161 (Quercus spp.) with respect to present and future climatic conditions. Molecular Ecology, 1162 25(23), 5907–5924. 1163 Rhoné, B., Defrance, D., Berthouly-Salazar, C., Mariac, C., Cubry, P., Couderc, M., Dequincey, 1164 A., Assoumanne, A., Kane, N. A., Sultan, B., Barnaud, A., & Vigouroux, Y. (2020). Pearl 1165 millet genomic vulnerability to climate change in West Africa highlights the need for regional 1166 collaboration. Nature Communications, 11(1), 5274. 1167 Robert, E., Lenz, P., Bergeron, Y., de Lafontaine, G., Bouriaud, O., Isabel, N., & Girardin, M. P. 1168 (2024). Future carbon sequestration potential in a widespread transcontinental boreal tree 1169 species: Standing genetic variation matters! Global Change Biology, 30(6), e17347. 1170 Robert, E., Lenz, P., Bergeron, Y., Isabel, N., & Girardin, M. P. (2025). Black spruce growth 1171 under climate extremes: Genetic insights for managing a key resource production species. 1172 Forest Ecology and Management, 597(123129), 123129. 1173 Rossi, S. (2015). Local adaptations and climate change: converging sensitivity of bud break in 1174 black spruce provenances. International Journal of Biometeorology, 59(7), 827–835. 1175 Rossi, S., & Bousquet, J. (2014). The bud break process and its variation among local 1176 populations of boreal black spruce. Frontiers in Plant Science, 5, 574. 1177 Schmidt, T. L., Jasper, M.-E., Weeks, A. R., & Hoffmann, A. A. (2021). Unbiased population 1178 heterozygosity estimates from genome/i1wide sequence data. Methods in Ecology and 1179 Evolution, 12(10), 1888–1898. 1180 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 47 Shafer, A. B. A., Cullingham, C. I., Côté, S. D., & Coltman, D. W. (2010). Of glaciers and 1181 refugia: a decade of study sheds new light on the phylogeography of northwestern North 1182 America. Molecular Ecology, 19(21), 4589–4621. 1183 Sniderhan, A. E., Mamet, S. D., & Baltzer, J. L. (2021). Non-uniform growth dynamics of a 1184 dominant boreal tree species (Picea mariana) in the face of rapid climate change. Canadian 1185 Journal of Forest Research, 51(4), 565–572. 1186 Urban, M. C., Bocedi, G., Hendry, A. P., Mihoub, J.-B., Pe’er, G., Singer, A., Bridle, J. R., 1187 Crozier, L. G., De Meester, L., Godsoe, W., Gonzalez, A., Hellmann, J. J., Holt, R. D., Huth, 1188 A., Johst, K., Krug, C. B., Leadley, P. W., Palmer, S. C. F., Pantel, J. H., … Travis, J. M. J. 1189 (2016). Improving the forecast for biodiversity under climate change. Science (New York, 1190 N.Y.), 353(6304), aad8466–aad8466. 1191 Vallejos, R., Osorio, F., & Bevilacqua, M. (2020). Spatial relationships between two 1192 georeferenced variables: With applications in R (2020th ed.). Springer Nature. 1193 Verrico, B. M., Capblancq, T., Fitzpatrick, M. C., & Keller, S. R. (2025). Reciprocal evaluation of 1194 genomic offset predictions of climate maladaptation with independent empirical datasets. 1195 The American Naturalist, 739111. https://doi.org/10.1086/739111 1196 Walker, X. J., Mack, M. C., & Johnstone, J. F. (2015). Stable carbon isotope analysis reveals 1197 widespread drought stress in boreal black spruce forests. Global Change Biology, 21(8), 1198 3102–3113. 1199 Zheng, X., Levine, D., Shen, J., Gogarten, S. M., Laurie, C., & Weir, B. S. (2012). A high-1200 performance computing toolset for relatedness and principal component analysis of SNP 1201 data. Bioinformatics (Oxford, England), 28(24), 3326–3328. 1202 1203 1204 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 48 Figure legends 1205 Figure 1. Sampling and population structure of black spruce. A. Sampling and location of 1206 common gardens. Dots show the location of populations colored by geographic regions 1207 delineated at equal intervals along the longitudinal axis. Black triangles indicate the location of 1208 common gardens (PR - Peace River, ML, Mont-Laurier, CH - Chibougamau, AC - Acadia). 1209 White area shows the distribution of black spruce. B. Pie charts showing mean genetic ancestry 1210 proportions per population with size of the piecharts proportional to population sample size. C. 1211 Admixture analysis showing ancestry proportions of samples at K=6 (including red spruce). D. 1212 Principal component analysis computed on 29,978 biallelic SNPs. Colors correspond to 1213 geographical regions depicted in A. E. FST matrix for five genetic clusters of black spruce. F. 1214 Observed heterozygosity and the number of polymorphic sites in black spruce populations. 1215 Populations are colored by geographical regions depicted in A. 1216 Figure 2. Local adaptation of black spruce to climate. A. Ordination of populations based on 43 1217 climate variables. Vectors of five climate variables identified with Gradient Forest are shown 1218 with arrows, with lengths scaled by a factor of 12. Locations of four common gardens (PR, ML, 1219 CH, AC) are shown with triangles. B. Importance ranking of top climate and spatial variables 1220 associated with genetic variation across black spruce distribution. C. Population means with 1221 standard deviation of height, biomass increment and survival rate plotted against climate 1222 transfer distance - Euclidean distance between the climate of origin of each planted population 1223 and that of the common garden site. D. Correlation coefficients between population means of 1224 traits and climate transfer distance. Bold font indicates coefficients statistically significant at P-1225 value<0.05 after multiple test correction. BI - biomass increment, Rc - recovery, Rr - relative 1226 resilience, Rs - resilience, Rt - resistance. Y ear following BI- labels indicates the beginning of 1227 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 49 the five-year period over which biomass increment was estimated. Colors in A and C 1228 correspond to geographical regions shown in Fig. 1A. 1229 Figure 3. Schematics of the validation procedure of genomic offset. A. Sampling scheme of 1230 populations used for training (circles) and testing (squares) in each common garden. When 1231 possible the same population was included in both train and test datasets, unless genotypes 1232 were not available (notably in AC). Triangles indicate location of common gardens. B. Steps 1233 required for model validation in a mock common garden example. Samples were partitioned into 1234 test and train datasets, then Gradient Forest models (GF) were trained on train datasets 1235 (population allele frequencies) and population-specific climate variables, and finally population 1236 genomic offsets (GO) were predicted for the common garden climate. For validation, genomic 1237 offset was compared against population means of fitness traits from test datasets to obtain 1238 Spearman rho, where more negative values indicated better predictive power of a given model. 1239 Numbers indicate investigation of the impact of six aspects of the data on model performance: 1240 1) impact of eight types of marker sets (variants) on model training, 2) comparison of 16 fitness 1241 traits with genomic offset, 3) impact of data limitations in train datasets on training and in test 1242 datasets on offset-fitness correlations, 4) comparing fitness-offset correlations for all populations 1243 (pops) vs. cluster-specific populations only, 5) impact of using populations from a single cluster 1244 vs. species-wide populations for training, 6) impact of using different sets of climate variables in 1245 training GF and predicting GO. 1246 Figure 4. Validation of genomic offset with height and biomass increments measured in 1247 common gardens. Spearman rho measures correlation between population genomic offset and 1248 population trait means, with more negative values indicating better correlation with fitness traits 1249 and therefore better model performance. A. Comparison of models trained on eight different 1250 marker sets validated using phenotypes in the train or test datasets. B. Comparison of models 1251 trained on eight different marker sets with genomic offsets validated against fitness separately 1252 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 50 within each genetic cluster. Dashed black lines indicate Spearman rho calculated between 1253 fitness trait and climate transfer distances. The marker types included subsets of 100 (0.1K), 1254 1000 (1K) or all available markers (~29K), subsets including low-frequency variants (0.1K-lf, 1K-1255 lf), and markers found to be in association with climate (LFMM, RDA, RDA-struct). 1256 Figure 5. Testing the impact of train dataset on model performance, where models were trained 1257 on complete datasets (“combined”) or datasets composed of populations from individual 1258 clusters. A. Models trained using three climate PCs as environmental variables. B. Models 1259 trained using selected five climate variables as environmental variables. 1260 Figure 6. Genomic offsets calculated with two Gradient Forest models for four 20-year periods 1261 and the medium emission scenario (SSP2-4.5). A. Results of genomic offset projected across 1262 four time periods using a model trained on the West cluster dataset and top five climate 1263 variables. B. Results of genomic offset projected across four time periods using a model trained 1264 on the combined dataset and three climate PCs. 1265 1266 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 51 Figures 1267 1268 1269 Figure 1. Sampling and population structure of black spruce. 1270 1271 1272 1273 51 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 52 1274 Figure 2. Local adaptation of black spruce to climate. 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 52 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 53 1285 Figure 3. Schematics of the validation procedure of genomic offset. 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 53 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 54 1296 Figure 4. Validation of genomic offset with height and biomass increments measured in 1297 common gardens. 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 54 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 55 1310 Figure 5. Testing the impact of train dataset on model performance, where models were trained1311 on complete datasets (“combined”) or datasets composed of populations from individual1312 clusters. 1313 1314 1315 1316 1317 1318 1319 1320 55 ed al .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 56 1321 Figure 6. Genomic offsets calculated with two Gradient Forest models for four 20- year periods1322 and the medium emission scenario (SSP2-4.5). 1323 1324 1325 1326 1327 1328 1329 56 ds .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 57 Supplementary materials 1330 Supplementary Figure 1. Patterns of isolation by distance across whole black spruce 1331 distribution. 1332 Supplementary Figure 2. Three major axes of spatial autocorrelation. 1333 Supplementary Figure 3. Spatial visualization of sPCA three main axes. 1334 Supplementary Figure 4. Observed heterozygosity versus number of polymorphic sites across 1335 populations of black spruce. 1336 Supplementary Figure 5. Performance of models trained on eight sets of markers validated 1337 against test and train datasets and all fitness traits. 1338 Supplementary Figure 6. Matrices of pairwise correlations between genomic offsets estimated 1339 using models trained on different sets of markers. 1340 Supplementary Figure 7. Matrices of pairwise correlations between genomic offsets estimated 1341 using models trained on different sets of markers. 1342 Supplementary Figure 8. Observed heterozygosity calculated for populations using selected 1343 marker sets. 1344 Supplementary Figure 9. Testing the impact of test sample size on model performance. 1345 Supplementary Figure 10. Testing the impact of removing populations from the train dataset on 1346 model performance. 1347 Supplementary Figure 11. Testing the impact of train set sample size on model performance. 1348 Supplementary Figure 12. Performance of models trained on eight sets of markers validated 1349 against three genetic clusters and all fitness traits. 1350 Supplementary Figure 13. Comparison of fitness population means with climate transfer 1351 distance across subset of fitness traits. 1352 Supplementary Figure 14. Comparison of fitness population means with climate transfer 1353 distance across subset of fitness traits. 1354 Supplementary Figure 15. Performance of models trained on single genetic clusters and three 1355 climate PCs in predicting fitness traits of the same or other genetic clusters. 1356 Supplementary Figure 16. Matrices of pairwise correlations between genomic offsets estimated 1357 using models trained on single genetic clusters or combined datasets and three climate PCs. 1358 Supplementary Figure 17. Performance of models trained on single genetic clusters and five 1359 selected climate variables in predicting fitness traits of the same or other genetic clusters. 1360 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint 58 Supplementary Figure 18. Matrices of pairwise correlations between genomic offsets estimated 1361 using models trained on single genetic clusters or combined datasets and five selected climate 1362 variables. 1363 Supplementary Figure 19. Matrices of pairwise correlations between genomic offsets estimated 1364 using models trained on two different climate datasets. 1365 Supplementary Figure 20. Change of climate across Canada in the current century. 1366 Supplementary Figure 21. Example of a confounding effect of differences in mean biomass 1367 increment between clusters on genomic offset estimates. 1368 1369 Supplementary Table 1. Geographical location of all initial black spruce populations 1370 (provenances) and the number of retained tree genotypes in each common garden after 1371 filtering. 1372 Supplementary Table 2. Geographical location and mean climate conditions of common 1373 gardens. 1374 Supplementary Table 3. List of climate variables. 1375 Supplementary Table 4. Number of populations assigned to test and train datasets. 1376 Supplementary Table 5. Number of genotypes assigned to test and train datasets. 1377 Supplementary Table 6. Results of AMOVA. 1378 Supplementary Table 7. Results of partial RDA analysis 1379 1380 Supplementary Dataset 1. Tree sample information. 1381 Supplementary Dataset 2. Sampling coordinates of red spruce populations. 1382 Supplementary Dataset 3. Past climate normals for 43 climate variables. 1383 Supplementary Dataset 4. Population means and standard errors of measured phenotypic traits. 1384 1385 1386 1387 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted October 31, 2025. ; https://doi.org/10.1101/2025.10.30.685617doi: bioRxiv preprint