Study of the Genetic Adaptation Mechanisms of Siberian Larch (Larix sibirica Ledeb.) Regarding Climatic Stresses Based on Dendrogenomic Analysis

Abstract

Dendrogenomics is a new interdisciplinary approach that allows joint analysis of dendrological and genomic data and opens up new ways to study the temporal dynamics of forest treelines, delineate spatial and temporal population structures, decipher individual tree responses to abiotic and biotic stresses, and evaluate the adaptive genetic potential of forest tree populations. These data are needed for the prediction of climate change effects and mitigation of the negative effects. We present here an association analysis of the variation of 27 individual tree traits, including adaptive dendrophenotypes reflecting the individual responses of trees to drought stress, such as the resistance (Rt), recovery (Rc), resilience (Rs), and relative resilience (RRs) indexes measured in 136 Siberian larch trees in 5 populations in the foothills of the Batenevsky Ridge (Kuznetsk Alatau, Republic of Khakassia, Russia), with variation of 9742 SNPs genotyped using ddRADseq in the same trees. The population structure of five closely located Siberian larch populations was relatively weak (F_ST = 0.018). We found that the level of individual heterozygosity positively correlated with the Rc and RR indices for the five studied drought periods and partly with the Rs indices for three drought periods. It seems that higher individual heterozygosity improves the adaptive capabilities of the tree. We also discovered a significant negative relationship between individual heterozygosity and the Rt index in four out of five periods, which means that growth slows down during droughts more in trees with higher individual heterozygosity and is likely associated with energy and internal resource reallocation toward more efficient water and energy usage and optimization of larch growth during drought years. We found 371 SNPs with potentially adaptive variations significantly associated with the variation of adaptive dendrophenotypes based on all three different methods of association analysis. Among them, 26 SNPs were located in genomic regions carrying functional genes: 21 in intergenic regions and 5 in gene-coding regions. Based on the obtained results, it can be assumed that these populations of Siberian larch have relatively high standing adaptive genetic variation and adaptive potential underlying the adaptations of larch to various climatic conditions.

1. Introduction

Droughts have become more frequent and drastic due to climate change [1,2], especially in boreal forests [3,4]. Therefore, it is very important to find whether tree populations possess sufficient standing genetic variation that may help trees to adapt to environmental stresses such as droughts. One of the most efficient approaches to find such variation is to search for genes whose variation is associated with the variation of adaptive dendrophenotypes using genome-wide genetic markers [5]. Our main hypothesis was that there are genes whose variation is significantly associated with individual responses of trees to the droughts, and therefore, our main objective was to study the association between the phenotypic variation of important adaptive traits related to the drought response and the genetic variation of individual trees and populations of Siberian larch (Larix sibirica Ledeb.) in the south area mostly affected by climate change [6,7,8,9,10,11].

Siberian larch is a coniferous tree common in Russia, Kazakhstan, Mongolia, and China. This is one of the key species of Siberian boreal forests, whose wood is superior in performance to many other coniferous species and has high ecological and economic significance. The wide range of Siberian larch includes both tundra zones in the north and forest and forest-steppe zones of Altai and Sayan Mountains in southern East Siberia [12,13,14].

The constant moisture deficit typical of the forest-steppe zone of southern East Siberia is gradually increasing with an increase in average annual temperatures due to global climate change, which, along with periodic droughts, affects the structure and dynamics of forest ecosystems. L. sibirica is a relatively drought-resistant species because of its high level of phenotypic plasticity and genetic variation. However, under the influence of rapidly increasing environmental stress, such as a decrease in the amount of summer precipitation, intense and recurring droughts, fires, and pest infestations, there is a concern that larch will not be able to effectively adapt to the new conditions associated with global warming [10,15,16].

Dendrochronological data based on measurements of annual tree wood growth rings contain a lot of information about the individual responses of trees to biotic and abiotic environmental factors [7,8,17,18,19,20,21,22,23]. In turn, the ability to individually genotype a large number of trees in populations for thousands of genes and genetic markers makes it possible to link genetic variation with the variation of specific dendrophenotypes [5,24,25,26,27,28]. This interdisciplinary approach, allowing the integration of dendrochronology, dendroecology, and genomics, has been gaining popularity in recent years [5,29,30,31].

According to the results of dendroclimatic studies conducted in the region of southern East Siberia, Siberian larch has successfully adapted to the modern conditions of the forest-steppe zone [20]. To understand the genetic mechanisms of this adaptation, we searched for DNA markers whose variation is associated with the variation of individual dendrophenotypes, reflecting the individual response of a tree to abiotic stress such as drought. For this purpose, genome-wide genotyping of a large number of trees in five populations was performed, and the structure of Siberian larch populations, genetic diversity, and genetic basis of adaptation to growing conditions within the forest-steppe zone of the foothills of the Batenevsky Ridge (Kuznetsk Alatau, southern East Siberia, Russia) were analyzed.

2. Materials and Methods

2.1. Sampling and Dendrophenotypes

The samples were collected from 150 trees in five populations in the foothills of the Batenevsky Ridge in the eastern part of the Kuznetsk Alatau, Republic of Khakassia (southern East Siberia, Russia) (Figure 1). This mountain system is characterized by a slight difference in elevation (500–1200 m above sea level), and most of its area is covered with mixed forest consisting of Siberian larch, Scots pine (Pinus sylvestris L.), and silver birch (Betula pendula Roth.). On the drier southern and southeastern slopes, open forest is interspersed with steppe areas. The forests of the study area are part of the state forest; forest districts are located on the ridge. The climate of the region is sharply continental, with large seasonal and daily temperature differences, hot summers, and frosty winters with little snow lasting from November to March. Despite the significant size of the study area, the high and similar variation in larch growth allows us to consider the forest-steppe ecotone of the Batenevsky Range as an integral geographical object, which is characterized by a unity of fluctuations in environmental conditions and the responses of vegetation to them [14]. Samples of larch wood and needles were collected at five population sites located 30–50 km apart in the forest-steppe zone along the foothills of the Batenevsky Ridge (Figure 1). The sites were selected on well-lit slopes oriented from southwest to southeast (

Table 1. The five population sampling sites of Siberian larch used in the study.

Population (Abbreviation Used in the Study)	Description	Coordinates	Altitude above Sea Level, m
Tuim (TUI)	Individuals and groups of larch trees along the steppe vegetation	54°24′20″ N 89°57′27″ E	550–600
Son (SON)	Mixed larch and birch forest	54°21′55″ N 90°22′04″ E	530–600
Bograd (BOG)	Mixed larch and birch forest	54°15′58″ N 90°41′30″ E	550–620
Bidja (BID)	Mixed larch, pine, and birch forest, with individuals and groups of larch trees in the steppe	54°00′20″ N 90°58′52″ E	640–670
Kamyziak (KAM)	Mixed larch and birch forest	53°55′52″ N 90°37′30″ E	700–770

).

Wood core samples were extracted from 150 mature, undamaged, live larch trees at chest level. The samples were collected and processed using standard dendrochronological techniques [32]. Individual tree ring widths (TRWs) were measured using the LINTAB device (Rinntech-Metriwerk GmbH & Co. KG, Heidelberg, Germany; https://rinntech.info/products/lintab, accessed on 14 October 2023) using the TSAP-Win program [33]. The TRWs were cross-dated, and the dating was checked using the COFECHA program [34]. The average TRW, generalized raw chronologies (avTRW) and their variance (varTRW) were used as dendrophenotypes.

Figure 1. Larch distribution area in Eurasia according to [14] (upper part, adapted with permission from Figure 1 in [13]) and location of the five population sampling sites of Siberian larch used in this study (lower part, adapted with permission from [35]): Bidja (BID), Bograd (BOGR), Kamyziak (KAM), Son (SON), and Tuim (TUIM). See also Table 1 for details.

Figure 1. Larch distribution area in Eurasia according to [14] (upper part, adapted with permission from Figure 1 in [13]) and location of the five population sampling sites of Siberian larch used in this study (lower part, adapted with permission from [35]): Bidja (BID), Bograd (BOGR), Kamyziak (KAM), Son (SON), and Tuim (TUIM). See also Table 1 for details.

Five drought periods (1951, 1963–1965, 1974–1976, 1999, and 2015, respectively) that occurred more than 6 years apart were selected from climatic extremes that significantly suppressed larch growth on a regional scale. To select drought periods, years with low values (<mean − SD) for the standardized precipitation–evapotranspiration index (SPEI) and/or precipitation averaged for the current April–July and/or previous June–September, as well as high temperatures (>mean + SD) over the current May–July and/or previous July-September, were identified in the study area (53.5–54.5° N, 89.5–91.5° E). Then, those periods with the simultaneous occurrence of several of these climatic deviations that also were accompanied by the formation of narrow tree rings (<mean − SD for averaged local tree-ring chronology) across all or most of the sampling sites were selected as drought periods. More detailed explanations of the selection of the drought periods exactly in these sampling sites, as well as the reasoning for the used seasonal timeframe of the climatic factors impacting significantly larch growth, can be found in [35]. To analyze the effect of this stress on larch growth, we used the resistance (Rt = Gd/Gprev), recovery (Rc = Gpost/Gd), and resilience (Rs = Gpost/Gprev) indices proposed by Lloret et al. (2011) [36], where Gprev is the average growth during the 3 years before the drought, Gd is the average growth during the drought, and Gpost is the average growth during the 3 years after the drought, as calculated from the tree-ring series of each tree after their standardization (omission of age-related trends fitted by cubic smoothing splines) using the ARSTAN program [37]. A more detailed description can also be found in [35]. The relative resilience indices (RRs) were also used, calculated as RRs = Rs − Rt. Additionally, the age of the tree at the time of collection (Age), the average length of 15 needles (avLn), and the variance of these 15 needle length measurements (varLn) were used as phenotypic traits. The radial growth trends (trendTRW) for each tree and the mean chronology for each population were calculated as the slope of the linear regression of the respective indexed time series after standardization (i.e., after elimination of the long-term trends related to tree size, age, and forest stand dynamics) for the 30-year period of the years 1990–2019. Phenotypic data that were not normally distributed according to the Shapiro–Wilk test were normalized using the rank-based inverse normal transformation (INT) algorithm. A total of 26 dendrophenotypic traits were considered in this study (

Table 2. Description of the dendrophenotypes.

Dendrophenotype	Abbreviation
Age of the tree at the time of collection	Age
Average needle length	avLn
Variance of needle length	varLn
Average tree ring width	avTRW
Variation of tree ring width	varTRW
Radial growth trends (1990–2019)	trendTRW
Index	Drought years
	1951	1963–1965	1974–1976	1999	2015
Resistance	Rt1	Rt2	Rt3	Rt4	Rt5
Recovery	Rc1	Rc2	Rc3	Rc4	Rc5
Resilience	Rs1	Rs2	Rs3	Rs4	Rs5
Relative resilience	RRs1	RRs2	RRs3	RRs4	RRs5

).

A pairwise Wilcoxon rank-sum test was performed to compare the dendrophenotypes and individual heterozygosity data between populations using the pairwise.wilcox.test command from the basic R package stats [38]. The Benjamini–Hochberg procedure was used as a method of multiple testing correction.

2.2. DNA Extraction and ddRAD Sequencing

DNA was isolated from the needles of the collected larch samples using the CTAB method [39]. The DNA concentration was assessed using a Qubit 2.0 fluorimeter (Thermo Fisher Scientific, Waltham, MA, USA). The purity and quality of the isolated DNA were also assessed using an Implen NanoPhotometer P330 (Implen, München, Germany). High-quality DNA samples with an A260/230 of approximately 1.8 and a concentration of 20–150 ng/μL were selected for this study.

The preparation of the ddRAD-seq libraries was based on a modified version of the protocol [40]. The DNA samples were treated with two restriction enzymes, EcoRI and MseI [41]. Validation of the prepared ddRAD-seq library pools before sequencing was performed on an Agilent 2200 TapeStation System (Agilent Technologies, Santa Clara, CA, USA). Single-end sequencing of the ddRAD-seq libraries was performed in three lanes using 100 cycles on a NovaSeq 6000 sequencer (Illumina, San Diego, CA, USA).

2.3. SNP Calling

The obtained sequencing data underwent several stages of initial processing: filtering, trimming according to quality indicators, and demultiplexing using the process_radtags utility included in the Stacks software version 2.65 [42]. Reads from each sample were then aligned to the Siberian larch reference genome [43] using Bowtie2 version 2.3 [44].

Variant calling was performed using the Gstacks utility from the Stacks software. The resulting set of alignment-covered loci was subjected to several filtering steps using the Populations utility to retain only loci that were present in at least 75% of all samples and in 75% of samples within each population. The maximum level of observed heterozygosity of the site should not exceed 0.85, and the frequency of the minor allele should not be less than 0.05. SNPs that were in linkage disequilibrium (LD) (r² > 0.8) were removed. Imputation of missing values was performed using the LD-kNNi method in TASSEL v.5.0 [45].

2.4. Genetic Structure of Populations

General indicators of genetic diversity, such as the mean observed (H_O) and expected (H_E) heterozygosities, nucleotide diversity (π), mean individual fixation index (F_IS), and their associated standard errors (s.e.), were estimated using the Stacks v2.5 software. The pairwise F_ST values among all the populations tested and their confidence intervals (using 10,000 bootstrap samples) were estimated using the R package StAMPP [46]. To determine the population structure, principal component analysis (PCA) was performed based on allele frequency data using the R ade4 package [47]. Sparse non-negative matrix factorization (sNMF) analysis was performed using LEA [48], and the estimation of individual ancestries was performed using the Admixture algorithm implemented in the AdmixPipe program [49]. Hierarchical analysis of molecular variance (AMOVA) was performed using Arlequin v. 3.5.1.2 [50]. The Mantel test to identify correlations between genetic and geographic distances was performed on two datasets: for 84 individual trees of the KAM, BID, and SON populations, for which the exact coordinates of the trees were recorded and the IBS genetic distance matrix was calculated; and for five populations using the average coordinates and the matrix of Nei’s genetic distances using the R package vegan v. 2.6-2 [51].

2.5. Associations between Dendrophenotypes and Individual Heterozygosity

For each sample tree, the individual level of heterozygosity was calculated as the number of polymorphic sites in the heterozygous state divided by the total number of polymorphic sites genotyped for a given sample. The correlations between individual heterozygosity and dendrophenotypes were analyzed using Spearman’s rank correlation and Pearson’s correlation coefficient using the R stats package [38].

2.6. Genotype–Phenotype Associations

The genotype–dendrophenotype associations were analyzed using the general linear model (GLM) and mixed linear model (MLM) in the TASSEL software version 5.0, where each site is an independent variable and each dendrophenotype is a dependent (response) variable. MLM analysis considers the confounding factors such as the pairwise kinship matrix K (reflecting relatedness between individuals) and population structure (the Q-matrix coefficients reflecting ancestry contributions into individuals), whereas GLM analysis accounts only for the population structure. The K-matrix of relatedness was calculated using the centered IBS method. The Q-matrix was obtained using the sNMF algorithm for the most probable number of population clusters (K), which equaled two for our samples. To search for associations whose weak signals were systematically observed in several different periods of drought, common overlapping SNPs were found for the SNP–dendrophenotype associations that passed the filtering threshold of p-value < 0.01.

The genotype–phenotype associations were also assessed using the Bayesian sparse linear mixed model (BSLMM) implemented in the GEMMA software package v.0.97 [52]. BSLMM is a polygenic model that accounts for the contribution of single SNPs with large effects, as well as the simultaneous contribution of multiple SNPs with smaller effects, on phenotypic variation. To achieve this, BLSMM includes the main effects of individual SNPs as predictors of the phenotype and the polygenic effect resulting from the combination of multiple SNPs with small effects. To identify significant genotype–phenotype associations for each trait, a posterior inclusion probability (PIP) threshold of >0.25 was used.

2.7. SNP Annotation

To analyze the genomic regions where SNPs associated with dendrophenotypes were located, annotation of the Siberian larch reference genome assembly was used [53]. Annotation of selected markers was performed using the SNPdat program [54]. Additionally, for SNPs associated with dendrophenotypes, the ±1000 bp long nucleotide sequences flanking the SNPs were extracted from the genome assembly using SNP coordinates, and a search for homologs of these extracted sequences was performed in the NCBI BLAST “nt” database [55].

3. Results

3.1. Dendrophenotypes

In order to understand how the values of dendrophenotypic traits differ between populations and among individuals, the pairwise Wilcoxon rank-sum test was performed (Table S1) and box-plots for each trait were generated (Figure 2), respectively. They demonstrated that some population samples were different from each other for some traits, and there was significant phenotypic variation between individuals (Table S1).

The average age (Age) of the studied trees was 124 years, ranging from 34 (kam2 tree) to 319 (bid24 tree) (Table S1). The highest average age was in the KAM sample, 153.8 ± 13.9 years; the lowest average age was in the SON sample, 102.8 ± 5.6 years (Table S1). Only between these two samples was there a significant difference in age based on the Wilcoxon rank-sum test (p = 0.019; Table S1).

The average needle length (avLn) for all the trees was 22.8 ± 0.3 cm, varying widely from 12.5 cm in the kam28 tree to 36.9 cm in the kam7 tree, and on average, from 20.1 ± 0.5 in the TUIM sample to 21.1 ± 0.5 in SON, 22.9 ± 0.5 in BID, 23.5 ± 0.5 in BOGR, and 26.5 ± 0.9 in KAM (Table S1). The samples differed significantly between each other (p < 0.05; Table S1), except for the BID–BOGR and SON–TUIM pairs.

The average variance of the needle length (varLn) across all the trees was 14.4 ± 0.8 cm (Table S1), varying insignificantly between samples from 10.5 ± 0.9 in BOGR to 12.2 ± 0.9 in TUIM, 12.7 ± 1.0 in SON, 14.3 ± 1.6 in BID, and 21.5 ± 3.0 in KAM, except only for the BOGR and KAM pair (p = 0.048; Table S1).

The average tree ring width (avTRW) for all the trees was 0.87 ± 0.08, varying widely among individual trees from 0.21 ± 0.12 in the kam28 tree to 2.53 ± 1.93 in the kam3 tree, and on average, in the samples from 0.77 ± 0.06 in TUIM to 0.78 ± 0.06 in BOGR, 0.83 ± 0.11 in KAM, 0.85 ± 0.10 in BID, and 1.10 ± 0.07 in SON (Table S1). Among all the sample pairs, only the SON sample was significantly different from all the other samples (p < 0.05; Table S1).

The variance of the tree ring width (varTRW) was varying also widely among individual trees, from 0.01 in the kam28 tree to 4.74 in the son21 tree (Table S1). On average, the largest varTRW was observed in the SON sample (0.97 ± 0.19), while in the other samples, it varied from 0.33 ± 0.04 in BOGR, 0.36 ± 0.07 in BID, 0.43 ± 0.06 in TUIM, and 0.63 ± 0.16 in KAM population samples (Table S1). The only SON sample was significantly different from the all the other samples (p < 0.05; Table S1).

The mean tree radial growth trend (trendTRW) values were negative in 99 out of 135 trees, varying also widely among individual trees, from −0.063 in the son4 tree to 0.040 in the bog9 tree (Table S1). On average, they were slightly negative in all the samples except BOGR (0.001 ± 0.003), varying from −0.022 ± 0.003 in SON to −0.012 ± 0.002 in TUIM, −0.008 ± 0.002 in KAM, and −0.004 ± 0.002 in BID (Table S1). The SON sample was significantly different from all the other samples, and the TUIM sample was also significantly different from the BOGR and BID samples (Table S1).

The growth practically completely ceased (Rt = 0) in 21 trees: bogr1, bogr3, bogr7, bogr9, bogr12, bogr14, kam17, and kam18 during both drought periods 1999 and 2015; tuim2 and tuim26 during the drought period of 1951; kam3, kam20, kam23, kam25, kam26, tuim8, tuim13, tuim15, tuim16, tuim26 during the drought period of 1999; and tuim14 during the drought period of 2015. The resistance (Rt) index very widely varied among individual trees, from Rt4 = 0.00173 in the bogr10 tree during the drought period of 1999 to Rt1 = 2.75 in the son30 tree during the drought period of 1951 (Table S1). Among all the drought periods studied, the highest average resistance index was observed for the drought period of 1963–1965 (Rt2 = 0.91 ± 0.03). Among all the population samples, the highest mean resistance index values were observed in the SON population sample for all the drought periods (Rt1 = 1.26 ± 0.10, Rt3 = 0.95 ± 0.04, Rt4 = 0.94 ± 0.06, Rt5 = 0.70 ± 0.07), except for the 1963–1965 period (Rt2 = 0.71 ± 0.03).

The recovery (Rc) index varied extremely widely among individual trees, from Rc1 = 0.26 in the son18 tree during the drought period of 1951 to as large a value as Rc4 = 290.04 in the bogr10 tree during the drought period of 1999 (Table S1). The highest average value of the Rc index was observed for the 1999 drought (8.31 ± 3.09). For the droughts of 1951 and 1974–1976, the highest Rc values were observed in the KAM population (Rc1 = 4.48 ± 1.31 and Rc3 = 1.40 ± 0.08, respectively); for the droughts of 1963–1965 and 1999, the highest Rc values were observed in the BOGR population (Rc2 = 1.33 ± 0.06 and Rc4 = 24.74 ± 15.38, respectively); and for the 2015 drought, the highest average value recovery index was observed in the TUIM population (Rc5 = 6.29 ± 1.34). The SON population showed the lowest Rc values (significantly different from the other populations) for all the drought periods.

The resilience (Rs) index varied very widely among individual, from Rs4 = 0.21 in the son23 tree during the drought period of 1999 to Rs5 = 2.83 in the kam26 tree during the drought period of 2015 (Table S1). The average value of the Rs index for 136 trees varied from Rs4 = 0.68 ± 0.02 for the 1999 drought to Rs5 = 1.33 ± 0.04 for the 2015 drought. Among the population samples, the highest average values of the Rs index were observed for trees of the KAM population for the droughts of 1951 and 1963–1965 (Rs1 = 1.01 ± 0.04 and Rs2 = 1.30 ± 0.07, respectively, Table S1) and for trees of the SON, BID, and TUIM populations for the droughts of 1974–1976, 1999, and 2015 (Rs3 = 0.89 ± 0.05, Rs4 = 0.85 ± 0.03, and 1.69 ± 0.07, respectively, Table S1). For three of the five drought periods, 1963–1965, 1999, and 2015, the lowest values of the Rs index were observed in the SON population (Rs2 = 0.60 ± 0.04, Rs4 = 0.48 ± 0.03, Rs5 = 0.90 ± 0.07, Table S1), which was significantly different from the other populations (Table S1).

The relative resilience (RRs) index varied also very widely among individual trees, from RRs1 = −1.75 in the son30 tree during the drought period of 1951 to RRs5 = 2.15 in the tuim1 tree during the drought period of 1951 (Table S1). The average value of the RRs index for all the trees varied from RRs5 = 0.78 ± 0.04 for the 2015 drought to RRs2 = 0.09 ± 0.02 for the 1963–1965 drought. The highest values of RRs were observed for the 1951 drought for trees of the KAM population (RRs1 = 0.59 ± 0.05), for the 1963–1965 drought for trees of the BOGR population (RRs2 = 0.25 ± 0.04), and for the last three droughts for trees of the TUIM population (RRs3 = 0.16 ± 0.03, RRs4 = 0.50 ± 0.04, and RRs5 = 1.20 ± 0.10, respectively, Table S1). The lowest values of the RRs index were observed in the SON population for all five drought periods (RRs1 = −0.50 ± 0.10, RRs2 = −0.11 ± 0.04, RRs3 = −0.06 ± 0.05, RRs4 = −0.46 ± 0.04, RRs5 = 0.20 ± 0.08, respectively, Table S1).

In general, the obtained dendrophenotypic data demonstrate that individual trees had very different individual responses to climatic stresses, such as drought, which well justifies the main objective of this study to associate the variation of individual dendrophenotypes that reflect individual tree responses to drought with the genetic variation of the large number of genetic markers genotyped in the same trees.

3.2. SNP Calling

For 150 larch tree samples, 1.8 billion 100 bp long reads were obtained via ddRADseq sequencing. Fourteen trees were excluded from further analysis because of insufficient sequencing quality. After the initial processing and quality filtering, approximately 1.4 billion reads with lengths from 32 to 92 bp were selected for further analysis, with an average of 10,496,058 ± 391,691 reads per sample. On average, 97% of the reads for each sample were successfully mapped to the Siberian larch reference genome. In total, 11,550 SNPs that met the filtering criteria were selected. For pairs of SNPs with LD (r² > 0.8), only one SNP was selected for further analysis. The final dataset contained 9742 biallelic SNPs genotyped in 136 Siberian larch trees (Data S1).

3.3. Genetic Structure of Populations

For each sample, general parameters of genetic variation, such as the nucleotide diversity (π), observed (H_O) and expected (H_E) heterozygosities, and fixation index (F_IS), were calculated and are presented in

Table 3. Values of mean nucleotide diversity (π), observed (H_O) and expected (H_E) heterozygosities, individual fixation index (F_IS), and their associated standard errors (s.e.) calculated for five population samples of Siberian larch based on all 9742 SNPs and 371 SNPs significantly associated with the variation of adaptive dendrophenotypes.

Population	π ± s.e.	HO ± s.e.	HE ± s.e.	FIS ± s.e.
9742 SNPs
TUI	0.00033 ± 0.00000	0.162 ± 0.001	0.162 ± 0.001	0.036 ± 0.017
SON	0.00030 ± 0.00000	0.118 ± 0.001	0.146 ± 0.001	0.163 ± 0.016
KAM	0.00033 ± 0.00000	0.147 ± 0.001	0.161 ± 0.001	0.099 ± 0.019
BOG	0.00034 ± 0.00000	0.173 ± 0.001	0.166 ± 0.001	0.012 ± 0.014 *
BID	0.00034 ± 0.00000	0.170 ± 0.001	0.170 ± 0.001	0.033 ± 0.015
All	0.00289 ± 0.00003	0.154 ± 0.001	0.166 ± 0.001	0.095 ± 0.049
371 SNPs
TUI	0.00189 ± 0.00012	0.171 ± 0.009	0.166 ± 0.007	0.029 ± 0.100 *
SON	0.00160 ± 0.00011	0.125 ± 0.008	0.141 ± 0.007	0.117 ± 0.020
KAM	0.00215 ± 0.00013	0.184 ± 0.009	0.189 ± 0.007	0.071 ± 0.101 *
BOG	0.00191 ± 0.00013	0.182 ± 0.010	0.168 ± 0.007	−0.010 ± 0.075 *
BID	0.00180 ± 0.00012	0.171 ± 0.009	0.159 ± 0.006	−0.005 ± 0.091 *
All	0.00192 ± 0.00011	0.166 ± 0.007	0.173 ± 0.005	0.073 ± 0.292 *

Note. All the values were significantly different from zero except the F_IS values marked by *.

. The levels of heterozygosity based on all 9742 SNPs was almost the same as those based on the 371 SNPs significantly associated with the variation of adaptive dendrophenotypes variation, although the latter ones were slightly higher if they were calculated for all the trees together as a single population.

The degree of genetic differentiation was assessed using the F_ST parameter (

Table 4. Pairwise F_ST fixation indexes with 95%confidence intervals in brackets calculated between 5 population samples.

Population	TUI	SON	KAM	BOG
SON	0.016 (0.014–0.017)
KAM	0.014 (0.012–0.016)	0.008 (0.007–0.009)
BOG	0.007 (0.005–0.008)	0.022 (0.019–0.024)	0.018 (0.015–0.020)
BID	0.007 (0.005–0.008)	0.020 (0.018–0.022)	0.017 (0.015–0.018)	0.007 (0.006–0.008)

). The maximum F_ST value was between the SON and BOG samples (0.022), and the minimum value between the TUI and BOG (0.007). The low pairwise F_ST values obtained indicate that the samples were relatively poorly differentiated and belong to genetically similar populations.

PCA was performed using allele frequency data to identify also the population structure. The analysis revealed that the studied samples do not form clearly defined unambiguous clusters, although two tentative clusters can presumably be distinguished: BID, BOG, and TUI form one cluster, and the SON and KAM samples form another (Figure 3).

In addition, for a better understanding of the population structure and further association analysis, a search for the most probable number of clusters or subpopulations K for K from 1 to 20 was conducted, with 20 replicates for each number K, and estimates of the contribution of each cluster to individual admixtures of trees (Q-values) were obtained. The ΔK method for selecting the most probable number K showed that the most probable numbers are K = 2 or K = 3. The contributions (“admixture”) of each cluster to individual trees (Q-values) according to the results of the Admixture program are shown in Figure 4. The sNMF algorithm, which was also used to calculate the Q-value matrix, showed similar results. When the number of clusters is K = 2, the trees belonging to the SON and KAM samples represent a separate cluster and are genetically different from the trees collected at other sites.

According to the Admixture analysis, all the genotypes were divided into two groups: the first group included the TUI, BOG, and BID populations, and the second group was SON–KAM.

AMOVA was conducted and Wright’s fixation indices (F-indices) were calculated for each hierarchical level (

Table 5. AMOVA results for five population samples of Siberian larch.

Source of Variation	Sum of Squares	Variance Components	Percentage Variation	F-Statistics over All Loci
Among groups	2079.3	7.7	1.0	FST = 0.018 FSC = 0.008 FCT = 0.010
Among populations within groups	3139.3	5.8	0.8
Within populations	195,233.5	731.7	98.2
Total	200,452.1	745.3	100.0

). The highest level of genetic variation was within populations (98.2%).

The correlation of the matrices of the pairwise geographic distances with the matrices of the pairwise genetic differences was calculated using the Mantel test: for 84 trees of the KAM, BID, and TUI samples and for 5 populations using average coordinates and Nei’s genetic distance matrix. The resulting test p-values indicate that there is no relationship between the geographic and genetic distances of populations (p = 0.171 and p = 0.558, respectively).

3.4. Associations between Dendrophenotypes and Individual Heterozygosity

For each sample, the individual level of heterozygosity was calculated, with an average 0.144 for all the samples. The lowest average values of individual heterozygosity were obtained for the SON population (0.112 ± 0.005) and the highest for the BOG population (0.164 ± 0.005). The Pearson’s and Spearman’s correlation coefficients of individual heterozygosity and each dendrophenotype were calculated (

Table 6. Pearson’s product-moment (r) and Spearman’s rank (r_S) correlation coefficients between dendrophenotypes and individual heterozygosity.

Dendrophenotype	r	p	rS	p
Age	−0.006	0.9413	0.011	0.8979
avLn	0.047	0.5858	0.030	0.7325
varLnr	−0.055	0.5262	−0.150	0.0808
avTRW	−0.049	0.5718	−0.046	0.5911
varTRW	−0.118	0.1701	−0.164	0.0560
trendTRW	0.372 *	<0.0001	0.371 *	<0.0001
Rc1	0.322 *	0.0004	0.366 *	<0.0001
Rc2	0.345 *	<0.0001	0.351 *	<0.0001
Rc3	0.339 *	<0.0001	0.437 *	<0.0001
Rc4	0.307 *	0.0007	0.435 *	<0.0001
Rc5	0.318 *	0.0003	0.331 *	0.0002
Rs1	0.064	0.4918	−0.028	0.7615
Rs2	0.228 *	0.0111	0.278 *	0.0019
Rs3	0.074	0.4098	0.062	0.4898
Rs4	0.180 *	0.0370	0.311 *	0.0002
Rs5	0.142	0.1097	0.200	0.0237
Rt1	−0.432 *	<0.0001	−0.378 *	<0.0001
Rt2	−0.076	0.4014	0.021	0.8144
Rt3	−0.286 *	0.0011	−0.258 *	0.0034
Rt4	−0.342 *	<0.0001	−0.280 *	0.0010
Rt5	−0.203 *	0.0216	−0.215 *	0.0146
RRs1	0.471 *	<0.0001	0.316 *	0.0005
RRs2	0.378 *	<0.0001	0.353 *	<0.0001
RRs3	0.405 *	<0.0001	0.493 *	<0.0001
RRs4	0.423 *	<0.0001	0.401 *	<0.0001
RRs5	0.216 *	0.0144	0.258 *	0.0033

* Significant values.

). These coefficients were calculated both for data not standardized by chronology and for indices calculated after the standardization procedure. The results were very similar; therefore, Table 6 presents coefficients only based on the unstandardized data.

There was a significant positive relationship between the recovery index (Rc), relative resilience index (RRs), and individual heterozygosity for all five drought periods (p-value < 0.05). The resistance index (Rt) showed a significant negative correlation with individual heterozygosity for four droughts, excluding the drought period of 1963–1965 (Rt2). The resilience index (Rs) correlated according to both the Pearson’s and Spearman’s correlation coefficients with individual heterozygosity based on the 1963–1965 and 1999 droughts, but only based on Spearman’s rank correlation coefficient for the 2015 drought data.

Dot-plot graphs of individual heterozygosity (IndHet) and dendrophenotypes (Rt, Rc, Rs, RRs, varTRW, and trendTRW), together with the Pearson’s (r) and Spearman’s rank (r_S) correlation coefficients, regression coefficients (R²) and trendlines based on linear, quadratic, cubic and local polynomial regression fitting, are presented for significant cases in Figure S1.

The mean radial growth and mean radial growth trend (trendTRW) are presented in

Table 7. The mean radial growth and mean radial growth trend (trendTRW) over the past 30 years during the period of 1990–2019 in five population sites of Siberian larch.

Population	BID	BGD	KAM	SON	TUIM
trendTRW	−0.002	0.002	−0.009	−0.017	−0.009
Mean radial growth	0.973	1.017	0.982	1.042	1.058

. The positive trend values mean that the growth rate was increasing over the past 30 years during the period of 1990–2019, taking into account the size–age trend, while negative values mean opposite—that it was decreasing. If the mean value of the growth index is more than one, then the trees have grown on average faster over the past 30 years than in their entire life; if less than one, then slower. The growth rate was slightly decreasing in four out of five populations, despite the fact that the trees have grown on average faster over the past 30 years than in their entire life in three out of five populations. A significant positive correlation between the radial growth trend (trendTRW) and individual heterozygosity was discovered. That is, the trees with higher heterozygosity demonstrate less of a decrease in the growth rate.

3.5. Genotype–Dendrophenotype Associations

According to the BSLMM, the proportion of variation in dendrophenotypes explained by the available genotypes (PVE) averaged 0.77, varying from the minimum value of 0.39 for the Rt5 index to the maximum value of 0.98 for avTRW. At the same time, the proportion of genetic variation explained by the variants with a relatively large effect (PGE) was 0.44 on average and varied from 0.03 in Rc4 to 0.96 in avTRW. For 21 out of 26 dendrophenotypes, the BSLMM method found associated SNPs with PIP values >25% (the estimated frequency by which an SNP has a sparse effect in the MCMC). In total, 2523 SNPs with presumably adaptive variation were selected using the BSLMM method, 453 of which were associated with more than 1 dendrophenotype.

According to the MLM, 1440 SNPs in total were associated with at least 1 of 26 dendrophenotypes, among which 565 were associated with 2 to 8 dendrophenotypes simultaneously.

According to the GLM, 2102 SNPs in total were associated with at least 1 of 26 dendrophenotypes, among which 1179 were associated with more than 1 dendrophenotype and 1 SNP (LS.4109781.252) was associated with 10 different dendrophenotypes.

A total of 4109 significant SNPs were detected using any of the 3 methods. Among them, 965 SNPs were detected by both GLM and MLM, and 371 SNPs by all 3 methods (Figure 5).

3.6. SNP Annotation

For the 371 SNPs found by all 3 association analysis methods, we analyzed the genomic regions in which these SNPs were localized (Table S2). According to the genomic annotation data of Siberian larch, only 26 SNPs were located in scaffolds with annotated genome features, whereas the rest were located in unannotated scaffolds. The majority of the successfully annotated adaptive markers were located in intergenic regions; only five were located in the coding regions of genes, and three of them were nonsynonymous and resulted in amino acid substitutions (

Table 8. Annotation of 26 potentially adaptive SNPs.

SNP	Dendrophenotype (GWAS Method)	Location *	Gene	Function	Reference
LS.3651056.309	Rt1 (GLM, MLM), avTRW (BSLMM)	Exon	Cytochrome P450 89A2	Oxidoreductase activity	[56,57]
LS.3905950.5973	avTRW, varTRW (GLM, MLM), Rs4 (BSLMM)	Exon	Xyloglucan endotransglucosylase/hydrolase protein B	Cell wall construction of growing tissues	[58]
LS.4033175.6261	Rs1 (GLM, MLM), RRs1 (BSLMM)	Exon	Cellulose synthase-like protein E6	Cellulose biosynthetic process	[59]
LS.4447596.15042	Age (GLM, MLM), Rc1 (BSLMM)	Exon	Uncharacterized protein LOC116133164	-	[53]
LS.4817075.929	varLn (GLM, MLM), Rc5 (GLM, BSLMM), RRs5 (GLM)	Exon	Protein Brevis radix-like 4	Modulator of root growth	[60]
LS.3903843.2180	Rc2, Rc3, RRs2, RRs3 (MLM, GLM) Rs2 (GLM), Rs3 (BSLMM)	78	RING-H2 finger protein ATL3	Ubiquitin ligase	[61,62]
LS.4185574.3908	varLn (GLM, MLM), Rs3 (BSLMM)	325	Alpha-mannosidase 2	Production of complex-type glycans	[63]
LS.7774.12431	RRs4, Rs2 (GLM, MLM), RRs5 (GLM, BSLMM) Rc5, Rt5 (GLM), Age (BSLMM)	3951	O-fucosyltransferase 19-like	Component of the gibberellin signaling pathway	[64]
LS.10638.27578	avTRW (GLM, MLM, BSLMM), Age, Rc4 (GLM)	5283	GPALPP motifs-containing protein 1 (lipopolysaccharide-specific response protein)	Lipopolysaccharide-specific response	[53]
LS.56935.13023	Rs4 (GLM, MLM, BSLMM)	5957	Probable aldo-keto reductase 4	Glyphosate degradation	[65]
LS.21791.50056	RRs1, Rc5 (GLM, MLM), Rs3 (BSLMM)	6066	Unknown	-	[53]
LS.61616.10689	Rc1, Rt1 (GLM, MLM), Rs5 (BSLMM)	6109	Hypothetical protein	-	[53]
LS.24961.61616	RRs5 (GLM, MLM), Rt4 (BSLMM)	6200	Unknown	-	[53]
LS.7196.6722	RRs5, Rs5 (GLM, MLM), RRs2 (BSLMM)	8888	Hypothetical heat shock protein	-	[53]
LS.105004.9682	avLn, Rt2 (GLM, MLM), Rc3 (GLM, BSLMM), RRs1, Rc1, Rt3 (GLM)	12,733	Phospholipase A1-Igamma3, chloroplastic	Catalyzes the hydrolysis of phosphatidylcholine	[66]
LS.7184.3035	Rs1 (GLM, MLM), RRs2 (BSLMM)	12,939	Unknown	-	[53]
LS.7252.104584	Rc3 (GLM, MLM, BSLMM), RRs1, Rc1, (MLM, GLM), Rs1 (GLM)	15,386	PLAT/LH2 domain-containing lipoxygenase family protein isoform 2	Defense mechanisms against pathogens	[67,68]
LS.14907.46110	Rt2 (GLM, MLM, BSLMM), Rs1 (GLM)	15,426	E3 ubiquitin ligase BIG BROTHER-related-like protein	It may limit the duration of organ growth and ultimately organ size by actively degrading critical growth stimulators	[69,70]
LS.15117.58377	RRs2, Rc2 (GLM, MLM), Age (GLM), Rs4 (BSLMM)	34,970	Cysteine and histidine-rich domain-containing protein RAR1 (CHRD1)	Gene-mediated disease resistance	[71]
LS.30352.56788	Rs4 (GLM, MLM), Rs3 (BSLMM)	54,635	Unknown	-	[53]
LS.16524.116919	Rs1 (GLM, MLM), Rs5 (BSLMM)	58,871	Uncharacterized FCP1 homology domain-containing protein C1271.03c	-	[53]
LS.12450.3129	Rs2, avTRW, varTRW (GLM, MLM), avLn (BSLMM)	64,166	RING-H2 finger protein	Ubiquitin ligase	[61,62]
LS.1920.81278	Rc1 (GLM, MLM), Rt1 (GLM), Rt3 (BSLMM)	65,622	U-box domain-containing protein 52-like	Functions as an E3 ubiquitin ligase.	[72]
LS.5783.87994	Rt1 (GLM, MLM), Rc1 (GLM), Rs3 (BSLMM)	69,130	ATP-dependent DNA helicase RECG-like	Recombination and DNA repair	[73]
LS.2549.42867	Rc3 (GLM, MLM), RRs2 (BSLMM)	69,752	Hypothetical protein GW17_00026793	-	[53]
LS.16630.144091	varLn (GLM, MLM), avLn (BSLMM)	83,423	Adenylosuccinate lyase	Synthesis of purine nucleotides	[74]

* Location in gene or distance to the nearest gene (bp).

).

For all the intervals of the nucleotide sequences (±1000 bp from the SNP or the maximum possible interval when the range of 2000 bp exceeded the scaffold size limits) containing SNPs associated with dendrophenotypes, homologs were searched in the NCBI BLAST database. Highly homologous sequences were found for 57 intervals containing SNPs, and most matches were found with genomic sequences and mostly unknown mRNAs and hypothetical proteins of white spruce (Picea glauca), Sitka spruce (Picea sitchensis), and Loblolly pine (Pinus taeda). Only for eight intervals were homologous sequences containing annotated genes or other structural units of the plant genome found in the NCBI GenBank database. For example, for the SNP locus LS.4643215.2651, which is associated with the dendrophenotype avLn according to GLM and MLM, as well as with the Rc3 index according to BSLMM, a homology was found with the mRNA encoding multiprotein-bridging factor 1c (MBF1, XM_028068828) of cowpea (Vigna unguiculata). It is a transcription cofactor whose molecular function is to form a bridge between transcription factors and the basal transcription machinery. They are involved in developmental processes and stress responses [75].

SNP LS.3970711.165 associated with the dendrophenotypes Rc4 and RRs4 according to GLM and MLM and with the RRs4 index according to BSLMM is located in a locus homologous to the mRNA encoding copper methylamine oxidase-like (XM_022162594) of oilseed sunflower (Helianthus annuus), which is involved in the mechanisms of protection against pathogens [76].

SNP LS.890128.210, associated with Rs3 according to BSLMM and Rs2 according to MLM and GLM is located in a locus homologous to the gene encoding the MYB6 transcription factor of Monterey pine (Pinus radiata, KY196189), which is involved in the regulation of the biosynthesis of secondary metabolites [77].

SNP LS.3903843.2180, which is associated with five dendrophenotypes (Rc2, Rc3, Rs2, RRs2 and RRs3) according to GLM, Rc2, Rc3, RRs2 and RRs3 according to MLM, and with Rs3 according to BSLMM, is located in a locus homologous to the mRNA of RING-H2 finger protein ATL3-like of wild soybean (G. soja, XM_028350422). There are also larch annotation data for this SNP (Table 8).

SNP LS.4382060.16010, which is associated with dendrophenotypes RRs1, Rs1, Rs3, Rs5, and Rt3 according to GLM, Rs3 according to MLM, and Rc1 according to BSLMM, is located in a locus homologous to the mitochondrial gene encoding ribosomal protein S4 of eastern white pine (Pinus strobus, MN9653 08), which is involved in the formation of plant disease resistance [78].

SNP LS.29920.18803, which is associated with avLn, Rs1, and Rt4 according to GLM and MLM, Rs2 and Rt2 according to GLM, and Rc1 according to BSLMM, is located in a locus homologous to the gene encoding polygalacturonase of breadseed poppy (Papaver somniferum, XM_026530000), which is responsible for the degradation of cell wall pectin and is involved in the processes of organ aging and the response to biotic stress in plants [79].

Homologues related to retrotransposons were identified for two sequences. Thus, the interval containing SNP LS.9095140.133, which is associated with the dendrophenotype Rs4 according to GLM and MLM and with avLn according to BSLMM, was aligned to the sequence encoding the RNA-directed DNA polymerase of cabbage (Brassica oleracea, XM_013761928). The sequence containing SNP LS.4042787.3127, which is associated with dendrophenotypes RRS1, RC1, RS2, RT1 and RT3 according to GLM, with RRS1, RC1 and RT1 according to MLM, and with RS5 according to BSLMM, was homologous to the sequence of retrotransposons of maritime pine (Pinus pinaster, DQ394069) and Monterey pine (AJ004945).

4. Discussion

The results of the presented genome-wide analysis of the structure and genetic variation of natural populations of Siberian larch in the foothills of the Batenevsky Range are generally consistent with previous conclusions about the relatively weak population structure of closely located Siberian larch populations: the value of F_ST = 0.018 observed in this study is very close to the value of F_ST = 0.017 based on a larger dataset of 25,143 SNPs obtained earlier for Siberian larch in the Altai-Sayan region [80].

The genetic differentiations between two groups identified on the basis of population structure analysis (F_CT), between populations within groups (F_SC), and between all the populations (F_ST) are consistent with the results previously obtained for larch and based on a larger dataset of SNPs [80,81]. The largest proportion of the genetic variance, 98%, was within the studied samples, and only 1% of the variance was due to differences between groups.

The results of the PCA showed that the samples were weakly genetically differentiated, which was expected because of the relatively short distances between sampling sites (about 50 km). Some trees from some populations demonstrated higher genetic similarity to trees from other populations according to the results of the cluster analysis (Figure 3). This is not surprising and can be expected, since the distance between populations is only tens of kilometers, there are no clear gaps between them such as treeless areas or forests that do not include larch, and there are no pronounced geographical barriers. Thus, some gene flow and genetic exchange can be expected between these populations, especially considering that larch is a wind-pollinated and wind-propagated species with small and light winged seeds. It is possible that the observed distribution reflects the predominance of southwestern and western winds in the basin of this region in combination with landscape restrictions. In particular, the SON and KAM populations are located approximately upwind of the remaining populations. Minimal genetic differentiation was observed between the TUI and BOG populations. There was no correlation between the geographic distance and the relatedness of trees in the study area.

The lack of a significant correlation between individual heterozygosity and the avTRW and varTRW dendrophenotypes is consistent with those previously obtained by Babushkina et al. [24] for two populations of Siberian larch based on the eight most polymorphic microsatellite loci (SSRs), although in the case of varTRW, it was negative according to Spearman’s rank correlation coefficient (r_S = −0.164) and almost significant (p = 0.056), which may partly indicate that increased heterozygosity stabilizes individual development, reducing the variance of wood growth. There was also no correlation between individual heterozygosity and age, average needle length, or variance in needle length.

According to recently obtained preliminary data for 124 Siberian pine (Pinus sibirica Du Tour) trees growing in the same region 150–250 km to the south of study area, there was no significant correlation between individual heterozygosity based on 13,914 biallelic SNPs and avTRW, although a significant negative correlation of the varTRW with individual heterozygosity was also found (r_S = −0.303) (data not published).

The presence of a significant positive correlation between individual heterozygosity and the recovery (Rc) and relative resilience (RRs) indices for all the studied drought periods and partly the resilience (Rs) indices for three drought periods suggests that differences in the drought resistance of individual Siberian larch trees are associated with the level of individual heterozygosity, the increase in which improves the adaptive capabilities of the organism.

However, we also discovered a significant negative relationship between the level of heterozygosity and the resistance index (Rt), which directly reflects the efficiency of larch growth during the drought years in four out of five periods. Perhaps one of the adaptive mechanisms of larch involves a sharp reduction in growth during periods of drought as part of a survival strategy, which is more pronounced in trees with a higher level of individual heterozygosity. Whether this is a consequence of increased homeostasis or heterosis and overdominance is difficult to determine and requires additional research.

It is interesting to note that the confidence intervals for the correlation trends are narrower in the upper-to-middle intervals (0.14–0.15) of the heterozygosity values, and he trendlines based on the quadratic, cubic and local polynomial regression fitting better approximate the correlation between individual heterozygosity and dendrophenotypes with higher regression coefficients (Figure S1). It can be a signature of the “optimal” heterozygosity reflecting the stable individual growth and response to the stress factors and providing the maximum adaptation of an individual or a population to the environment in which the individual or population was formed and in which it now exists (in sensu [82,83,84]).

The search for genotype–phenotype associations using three different methods revealed 371 potentially adaptive SNPs. Only 26 SNPs were located in genomic regions carrying functional genes: 21 in intergenic regions and 5 in gene-coding regions. The genetic diversity summary parameters based on 371 potentially adaptive SNPs were only slightly higher than those that were based on all 9742 SNPs (Table 3). It seems that a limited number of SNPs are needed to appropriately assess genetic diversity.

It is important to note that the growth rate was slightly decreasing in four out of five populations over the past 30 years during the period of 1990–2019. This is likely related to the climate change, particularly reflecting the increased water deficit caused by warming. These observations are corroborated by recent climate-driven decreasing trends in the grain crop yields observed in the steppe territories of the study region [85]. However, creation of large water reservoirs has so far mitigated most of the drying climatic trends in the region [86] comparing to most of continental Asia, where severely declining tree growth and forest die-off were observed recently in many places on southern and lower fringes of forested areas [6,87].

5. Conclusions

Several regions of the genome and genes have been identified, the variation of which is associated with the variation of important adaptive traits, dendrophenotypes. We certainly realize the limits of the RADseq method and that the obtained genomic data could be insufficient to uncover all the genes and genetic mechanism that are responsible for adaptation to drought, although those that are found could be of interest for the scientific community. We identified several SNPs in candidate genes whose variation was associated with important adaptive dendrophenotypes. Based on these results, it can be assumed that Siberian larch has relatively high adaptive genetic variation and potential in the studied area. It seems well adapted to abiotic stresses such as droughts. However, the negative growth trend observed in four out of five populations could be a sign of the beginning negative effect of climate change. The results of this study will allow for a deeper understanding of the genetic mechanisms underlying the adaptations of larch to various climatic conditions. This study allowed us to detect important genes, and we obtained SNPs with significant adaptive variation that can be potentially used as genetic markers in the marker-aided breeding for stress tolerance as well as for creating a SNP genotyping chip for monitoring adaptive genetic variation in other larch populations, although more SNP genotyping data, preferably based on the whole-genome sequencing, are required for developing a comprehensive high-density SNP genotyping assay to study larch. The presented approach, which includes both dendrophenotypic and genetic data, can serve as a scientific basis for optimizing nature management, developing methods for the rational use of the studied species, identifying populations with good genetic potential and conducting environmental monitoring.