Population Structure, Genetic Diversity and Candidate Genes for the Adaptation to Environmental Stress in Picea koraiensis

Picea koraiensis is major silvicultural and timber species in northeast China, and its distribution area is an important transition zone for genus spruce migration. The degree of intraspecific differentiation of P. koraiensis is high, but population structure and differentiation mechanisms are not clear. In this study, 523,761 single nucleotide polymorphisms (SNPs) were identified in 113 individuals from 9 populations of P. koraiensis by genotyping-by-sequencing (GBS). Population genomic analysis showed that P. koraiensis was divided into three geoclimatic regions: Great Khingan Mountains climatic region, Lesser Khingan Mountains climatic region, and Changbai Mountain climatic region. Mengkeshan (MKS) population on the northern edge of the distribution area and Wuyiling (WYL) population located in the mining area are two highly differentiated groups. Selective sweep analysis showed that MKS and WYL populations had 645 and 1126 selected genes, respectively. Genes selected in the MKS population were associated with flowering and photomorphogenesis, cellular response to water deficit, and glycerophospholipid metabolism; genes selected in the WYL population were associated with metal ion transport, biosynthesis of macromolecules, and DNA repair. Climatic factors and heavy metal stress drives divergence in MKS and WYL populations, respectively. Our findings provide insights into adaptive divergence mechanisms in Picea and will contribute to molecular breeding studies.


Introduction
Picea is the third-largest genus in Pinaceae and a major component of coniferous forests in the northern hemisphere. Most conifers have the characteristics of long generations, large population size [1][2][3], and a low level of interspecific reproductive isolation [4]. Spruce underwent rapid differentiation during the Pliocene, Qinghai-Tibet Plateau uplift, and Quaternary glacial oscillation [5]. The continuous gene flow among sympatric and parapatric species makes the evolutionary relationship more complex [6,7]. Geographic isolation and gene introgression are two possible driving factors of species differentiation [8,9]. In conifers with long generations, a long period of time can elapse from the onset of differentiation to the cessation of gene flow. Studies have suggested that gene flow between different species continues to decline over time when the species is under either isolation or selective pressure [10]. Understanding gene introgression and selection pressure is of great significance for further understanding spruce phylogeny.
Picea koraiensis, belonging to the genus Picea, is distributed in the northeast of the Great Khingan Mountains, Lesser Khingan Mountains, Zhang Guangcai Mountains, Wanda (GF) and Zhanhe (ZH) populations. Group II and Group III are sister clades. The results of the phylogenetic tree showed that the genetic relationship of different populations of P. koraiensis was related to geographical distribution. The population genetic structure was further analyzed by ADMIXTURE v1.23 software [23]. The number of clusters (K) was set from 2 to 10, and the best number of clusters was selected according to the minimum cross-validation rate ( Figure 2). The optimal number of clusters was three. When K = 2 to 4, the error rate of cross-validation was small, indicating that 2 to 4 clusters is more reasonable. When K = 2, 127 samples were divided into P. koraiensis and the outgroup, and the MKS population of P. koraiensis had a small amount of outgroup gene infiltration. When K = 3, P. koraiensis was divided into 2 groups, the WYL population and other populations. When K = 4, other populations had extensive gene flow except for the WYL and GF populations. The results showed that gene introgression occurred in the MKS population and outgroups, and the WYL population had greater differentiation from other populations. To further analyze the population genetic structure, 112 individuals of P. koraiensis were analyzed by principal component analysis (PCA) (Figure 3). The results showed that P. koraiensis could be divided into three main clusters: MKS, WYL, and other populations. The results of PCA and admixture analyses indicated that there is an obvious differentiation between the two populations (MKS and WYL) and others.

Population Diversity and Selective Sweep Analysis
The data are shown in Table 1; the π(nucleotide diversity) of the P. koraiensis population was 0.00194 to 0.00231. Among the populations, the π value of the WYL population was the lowest (0.001943), while the π value of LJ was the highest (0.00235). The F st (fixation index) among P. koraiensis was 0.0254 to 0.1569, among which the mean F st between the WYL population and other populations was the highest (0.1211), followed by the MKS population (0.0983), and the mean F st between the TQL population and other populations was the lowest (0.0538). This result is consistent with the phylogenetic analysis, indicating that the WYL and MKS populations have a high degree of differentiation from other populations, and these two populations may be isolated from other populations and have less gene flow, resulting in a significant decrease in nucleotide diversity. That is, these two populations experienced a high degree of selection in the process of evolution. The degree of differentiation of the LJ, TQL, and ZH populations was low, and the three populations had frequent gene flow with other populations and high nucleotide diversity. The Tajima's D values of the nine populations of P. koraiensis were all positive (mean 0.1554), indicating that these nine populations may have experienced balanced selection.

Population Diversity and Selective Sweep Analysis
The data are shown in Table 1; the π(nucleotide diversity) of the P. koraiensis population was 0.00194 to 0.00231. Among the populations, the π value of the WYL population was the lowest (0.001943), while the π value of LJ was the highest (0.00235).

Population Diversity and Selective Sweep Analysis
The data are shown in Table 1; the π(nucleotide diversity) of the P. koraiensis population was 0.00194 to 0.00231. Among the populations, the π value of the WYL population was the lowest (0.001943), while the π value of LJ was the highest (0.00235). The main purpose of selective sweep analysis was to find the selected genes of the population, which further reveals the adaptive mechanism of population evolution. In our Plants 2023, 12, 1266 6 of 18 study, four populations of Changbai Mountain group were used as background group. The MKS and WYL populations were used as test groups, respectively. Fixation index (F st ) and nucleotide diversity ratios (π ratio) were used in combination to identify regions where selective sweeps occurred during species differentiation. A total of 645 and 1126 selected genes were identified in the MKS and WYL populations, respectively (Tables S1 and S2).

Geographical Isolation, Climate Heterogeneity, and Gene Introgression
In our study, whole-genome SNPs were obtained by the GBS-seq technique, and the phylogeny, population structure, and genetic differentiation of P. koraiensis distributed in northeast China were analyzed. The results showed that the genetic relationship among different populations of P. koraiensis was consistent with the division of geographical and climatic regions: there was gene introgression between the MKS population and outgroup, and there was obvious genetic differentiation among the MKS and WYL populations and other populations. Climate heterogeneity, gene introgression, and geographical isolation are the main reasons for the differentiation of P. koraiensis.
Climatic factors such as temperature and precipitation are likely to be among the important factors driving the differentiation of P. koraiensis due to different habitats. The phylogenetic results divided P. koraiensis into three groups corresponding to three geoclimatic regions: the Great Khingan Mountains climatic region, the Lesser Khingan Mountains climatic region, and the Changbai Mountain climatic region. The climate of the northern Great Khingan Mountains area is characterized by long sunshine, low precipitation, and low temperature; the climate of the Changbai Mountain region is characterized by medium sunshine, medium precipitation, and medium temperature; and the climate of the Lesser Khingan Mountains area is characterized by short sunshine, low precipitation, and low temperature [20]. Zhang et al. [20]  are also examples of temperature and moisture influencing population differentiation, such as Liaodong oak [24], eucalyptus in Australia [25], and alpine conifers in the western Alps [26].
Gene flow and introgression are also important causes of intraspecific differentiation in P. koraiensis. On the one hand, the size of gene flow between different groups within the species and the degree of differentiation were inversely proportional. These two groups, WYL and MKS, have little gene exchange with other groups, low nucleotide diversity, and greater differentiation; the ZH, TQL, and LJ populations have extensive gene flow with other groups, high nucleotide diversity, and low differentiation. On the other hand, gene introgression between species promotes the differentiation of different populations within P. koraiensis. The results of phylogenetic and population differentiation analysis showed that the MKS population was more differentiated from other populations and had a small genetic component of outgroups. Because of the large overlap in distribution between the MKS population of P. koraiensis and its sympatric species P. jezoensis, we speculate that there is a small amount of gene introgression from P. jezoensis to the MKS population. Many closely related species of P. koraiensis have been successfully crossed with closely related species of P. jezoensis in hybridization experiments [15]. Fossil evidence suggests that during the Pleistocene, P. jezoensis and its extinct relatives were widespread in northeastern Russia [27], that the ancestors of P. koraiensis and its relatives may have appeared in east Asia during the same period [28], and that the ancestors of P. jezoensis and P. koraiensis underwent infiltration before P. koraiensis and its relatives' migration to Eurasia [29]. Gene introgression between other species of the genus spruce has been reported repeatedly: P. glauca, P. engelmannii, and P. sitchensis in western North America [30][31][32][33][34][35] and P. likiangensis, P. purpurea, and P. wilsonii in the Qinghai-Tibet Plateau [6,7,[36][37][38] Geographical isolation is an important reason for the intraspecific differentiation of P. koraiensis. At the northern edge of the P. koraiensis distribution area, the MKS population has the characteristics of geographical and genetic isolation because of spatial isolation and limited gene communication; that is, the differentiation characteristics of the MKS population accord with the central-marginal hypothesis, in which marginal populations show lower genetic variation and higher differentiation than central populations. There are similar examples in other species, such as Abies sachalinensis in Hokkaido and Pinus cembra L. in the Alps [39,40]. The WYL population is distributed in the Lesser Khingan Mountains gold mining area. Due to the influence of mining and different geological conditions, the natural environment shows fragmentation and islanding, forming an ecologically isolated island, which is different from other distribution areas of P. koraiensis. Furthermore, there is geographical isolation. There are similar examples of other species distributed in mining areas, such as Caryophyllaceae [41,42], Scopelophila cataractae [43], and Colocasia esculenta [44]. In summary, geographical isolation leads to the differentiation of species and the formation of obvious geographical differences to a great extent; the same is true of Lamiaceae and Circaeasteraceae [45,46].

Molecular Mechanisms of Population Differentiation
To further detect natural selection signals in different populations of P. koraiensis and explore the differentiation mechanism, we performed selection sweep analysis to identify candidate genes based on KEGG and GO enrichment analyses involving various biological functions, using populations of the Changbai Mountain group as background populations and the more differentiated MKS and WYL populations as control populations. The key functional genes for which differentiation between different populations of P. koraiensis was subject to selection were to improve the resistance of the population, to promote flowering and fruiting ability for adaptation to different habitats, and to continuously reproduce offspring.
The MKS population is distributed in the northern part of the Great Khingan Mountains and is geographically distant from other populations, with large climatic differences. Climatic heterogeneity drives significant differentiation of the MKS Mountain population, and the function of genes shows the molecular adaptation mechanism of the MKS population to its habitat. The northern Great Khingan Mountains have long sunshine hours (approximately 2600 h per year and 2300 h per year in the Changbai Mountain region), and light affects the photoperiod, circadian rhythm, and photomorphic building of plants. The results of our selection sweep analysis showed that 15 genes in the MKS population were enriched in response to radiation, and 14 genes were enriched in response to light stimulus. Among these genes, CUL4 in Arabidopsis can regulate the flowering time in a short period of time during the day and is involved in photomorphogenesis [47,48]; GATA9 is involved in the regulation of circadian rhythms in Arabidopsis [49]; SBH2 is involved in photomorphogenesis [50]; and MBD9 regulates photoperiod and flowering in Arabidopsis through DNA methylation [51,52]. The northern Great Khingan Mountains are characterized by low annual precipitation (463 mm average annual rainfall and 680 mm average annual rainfall in Changbai Mountain), and our results show that 5 of the genes are enriched in response to water deprivation, with RPK2 being involved in pollen germination, response to cold, and response to water deprivation [53]; CYT1 being involved in defense responses to bacteria, hormone homeostasis, and salt stress [54,55]; and MGL being involved in cellular response to sulfate starvation and cellular response to water deprivation [56]. The northern Great Khingan Mountains are characterized by low temperatures (−2.4 • C annual mean temperature and 3.2 • C annual mean temperature in Changbai Mountain), and organisms living in this environment face various growth-related challenges from low temperatures, such as reduced lipid membrane fluidity. The lipid bilayer transmits external signals to its interior and protects the integrity of various organelle membranes in unfavorable environments. The results showed that 86 genes were associated with membrane-bound organelles and that 6 genes were involved in glycerophospholipid metabolism. Among these genes, PLC is involved in leaf membrane lipid remodeling [57] and regulation of stomatal movement and root development under abiotic stress [58,59]. Remodeling of membrane lipids may protect membranes from low temperatures. We hypothesize that these genes may play a major role in maintaining membrane integrity under low-temperature conditions. The functions of these genes suggest that MKS populations have evolved a range of genes adapted to the long light hours, low annual rainfall, and low temperatures of the habitat.
Among the Lesser Khingan Mountains group, the WYL population is distributed in the gold mining area, and there are also large amounts of iron ore and agate in its distribution area. The WYL population is distributed at a low altitude (average elevation approximately 350 m), and due to tailings pollution caused by mining, monsoon effects, and rainfall leaching, many heavy metals, such as Cd, Cr, Cu, Ni, Pb, and Zn, in the soil exceed background values, especially Pb, Cd, Zn, and Cu [60,61]. Heavy metals cause oxidative damage, disruption of cellular homeostasis, DNA strand breaks, fragmentation of proteins or cell membranes, and damage to light and pigments [62]. Plants undergo a series of physiological changes in response to heavy metal stress to adapt to their environment. The first line of defense of plants against heavy metal stress is to block the transport of metal ions in the roots through the plasma membrane [63]. The enrichment analysis showed that 19 genes are involved in the transport of metal ions, of which 5 genes were enriched in ABC transporter, where the ABC transporter family (e.g., ABCB1, ABCA15, ABC1, ABCG28, and ABC1K8) is involved in the efflux of metal ions from the plasma membrane and in the resistance of plants to heavy metals and aluminum toxicity [64,65]. For example, the ABCB family genes enhance resistance to cadmium and lead in Arabidopsis [65], and ABCB1 is upregulated in mango when it is subjected to lead stress [66]. After heavy metals penetrate the first line of defense and enter the cell, plants tolerate or neutralize the toxicity produced by heavy metals through the synthesis of biomolecules such as metal chelators and specific amino acids [63]. Forty-eight genes were identified as regulators of cellular macromolecular biosynthetic processes, among which PCYT2, PCYT1, EPT1, and PHO1 are involved in the synthesis of phosphate and hypophosphate [67][68][69][70][71][72][73]; ANS, FLS, and ANR are involved in flavonoid biosynthesis [74][75][76]; CysCK is involved in the synthesis of cysteine [77]; and G6PD, SPE3, GST, and RRM1 are involved in glutathione synthesis [78][79][80][81]. Phosphates, flavonoids, cysteines, glutathione, and other phytochelators are markers of heavy metal stress in plants [82,83]. GST genes are upregulated in horse rushes, tomatoes, and radishes when they are subjected to heavy metal stresses such as cadmium [84][85][86]. When both of the above mechanisms against heavy metal stress are exhausted, plants activate oxidative stress defense mechanisms and the synthesis of stress-related proteins and signaling molecules [63]. Seven genes were associated with oxidative stress defense in plants, among which MPK6 and FLS2 were enriched in the MAPK pathway and IAA, SAUR, and GID1 were enriched in the hormone signaling pathway. In Arabidopsis, MPK6 can be induced by Cd ions and Cu ions [87][88][89][90], and the IAA (AUX) concentration increases to promote lateral root growth and protect plants from Cd [91]. In addition, the regulation of transcription factor family expression by heavy metal stress has been reported [92][93][94]. The transcription factors identified in this study are members of the MYB family (MYB4, MYB52, MYB58, MYB78), WRKY family (WRKY7, WRKY35, and WRKY42), NAC, and BZIP68. In Arabidopsis, MYB family genes such as MYB4 are highly induced by Cd and Zn metal stress [95]. The expression of WRKY family members was significantly elevated under Cu and Cd metal exposure [96]. BZIP transcription factor expression was also induced by Cd stress [97]. For DNA damage caused by heavy metals, plants initiate DNA damage repair mechanisms with six genes involved in DNA repair, of which RBX1 and RFA1 are involved in nucleotide excision repair [98,99], MLH1 and PMS2 are involved in mismatch repair [69,100,101], and TOP3 and BRCA2 are involved in homologous recombination [102,103], suggesting that the WYL population evolved multiple DNA repair pathways. In addition, compared with Changbai Mountain, the sunshine in the Lesser Khingan Mountains area is short and the temperature is low. The WYL population also identified genes related to light morphogenesis, circadian rhythm, and cold tolerance. For example, in Arabidopsis thaliana, the CSN4 and CSN8 participate in the development of seedlings in the dark and promote the development of light morphology in the dark [104,105]; ABC1K8 and BLH1 are involved in regulating circadian rhythms [106,107]; and LSM1A, CIPK9, and PSRP2 are involved in plant adaptation to cold environment [108][109][110]. Our results show that the genome of the WYL population not only responds to the climatic environment of the Lesser Khingan Mountains area, but also may be stressed by heavy metals in mining areas, and the genome has undergone adaptive evolution.

Sample Collection and GBS Analysis
We collected samples according to the 4 climatic zones of P. koraiensis and identified a total of 9 populations of 10 to 15 individuals each, for a total of 112 individuals (Figure 7, Table 2 . Three species, P. jezoensis (sympatric species of P. koraiensis), P. likiangensis, and P. pungens, were used as outgroups, with five individuals per species. Plants from each population were separated by at least 100 m at the time of sampling, and needles were dried directly using silica gel and kept at low temperature. DNA extraction from plant needles was performed using a modified bromocetyltrimethylamine (CTAB) method, followed by a Nanodrop-1000 spectrophotometer (Nanodrop, MA, USA) to detect nucleic acid concentration and DNA purity on 1% agarose gel electrophoresis.

Sample Collection and GBS Analysis
We collected samples according to the 4 climatic zones of P. koraiensis and identified a total of 9 populations of 10 to 15 individuals each, for a total of 112 individuals ( Figure  7, Table 2 . Three species, P. jezoensis (sympatric species of P. koraiensis), P. likiangensis, and P. pungens, were used as outgroups, with five individuals per species. Plants from each population were separated by at least 100 m at the time of sampling, and needles were dried directly using silica gel and kept at low temperature. DNA extraction from plant needles was performed using a modified bromocetyltrimethylamine (CTAB) method, followed by a Nanodrop-1000 spectrophotometer (Nanodrop, MA, USA) to detect nucleic acid concentration and DNA purity on 1% agarose gel electrophoresis.    DNA samples were qualified, sequenced, and constructed in GBS libraries according to the methods reported by Poland et al. [111]. Each sample extracted from 1.5 µg of DNA was digested in 96-well plates with EcoRI and NiaIII restriction enzymes. The digestion reaction was performed in 96-well plates using 100 ng of DNA for each individual reaction digested with EcoRI and NIaIII (New England Biolabs, Ipswich, MA, USA). The digested product was mixed with A1 and A2 linkers at 25 pmoL per well, adding linkers at both ends of the DNA. Library pooling, selection of size (400-600 bp) on a 1% agarose gel, column purification using PCR purification kit (NEB), and amplification of 12 cycles using Phusion DNA polymerase (NEB) were carried out. After the second column wash, the average fragment size was estimated on a BioAnalyzer 2100 (Agilent, Santa Clara, CA, USA) using a DNA1000 chip, and library quantification was performed using PicoGreen (Invitrogen, Carlsbad, CA, USA). Sequencing with PE125 on HiSeq4000 (Illumina, San Diego, CA, USA) adjusted the common library to 10 nmol. Numerous raw sequencing sequences (raw reads) were obtained through high-throughput sequencing. To ensure data quality, the quality control of the original data should be carried out before information analysis, and data filtering should be used to reduce data interference information. We filtered the dismounted raw reads more strictly to obtain high-quality clean reads for subsequent information analysis. The conditions for filtering were as follows: (1) reads containing adapter sequences were removed; (2) reads containing unknown base N ratios greater than 10 were removed; and (3) low-quality reads (the number of bases of mass value Q ≤ 10 accounted for more than 50% of the entire read) were removed for genomic electronic digestion assessment. To check the efficiency of digestion and analyze all markers, we digested the reference genome according to the site of restriction enzymes and counted the corresponding conditions for later evaluation. Sequencing evaluation: statistics are made of output data, including sequencing data yield, sequencing error rate, Q20 content, Q30 content, N content, GC content, etc., to evaluate whether the library construction and sequencing are successful.
Using the comparison software Burrows-Wheeler Aligner (BWA) ver.0.7.12 (Li, Cambridge, U.K.) [112], the filtered reads were compared with the reference genome (European spruce) by the mem algorithm, and the alignment parameter was-k 32-M. After the comparison, the results were marked by Picard (1.129) software. The Unified Genotyper module of the software GATK (3.446) is used to detect the variant of multiple samples of the processed comparison files, and the detected variations are filtered by VariantFiltration. The filtering parameters are -Window 4, -filter "QD < 4.0 || FS > 60.0 || MQ < 40.0", -G_filter "GQ < 20". ANNOVAR was used to make functional comments on the detected variant. Using VCFTOOL v.0.1.11 (Danecek, Cambridge, UK) [113], the SNP loci (127 individuals) obtained by GBS sequencing were filtered if the missing ratio ≤20% dint maf ≥0.01 and all indels were removed.

Phylogenetic and Population Structure Analyses
The phylogenetic tree of 127 individuals was constructed by the Neighbor-Joining method and the P-diatance model [114] based on MEGA version X software [115] and corrected by Jukes-Cantor, Tajima-Nei, and the Maximum Composite Likelihood model. A bootstrap consensus tree was formed from a sample of 1000 replicates. iTOL version 6.4 was used to edit and beautify the phylogenetic tree [116].
To carry out the cluster analysis of the population structure, the parameters of the ancestry coefficient matrix Q and the population allele frequency matrix F and K were calculated by admixture software [23] and finally plotted by the R software package (http://www.r-project.org/, accessed on 12 January 2023). The K value range was set from 1 to 100, and the default settings were used for the other parameters. The optimal clustering number was determined according to the minimum value of the cross-validation (CV) error, and the R language package was used to draw the results.
The genetic structure of the spruce population was evaluated by principal component analysis (PCA). We used GCTA software (http://cnsgenomics.com/software/gcta/pca. html, accessed on 12 January 2023) to reduce the dimension of the selected SNP loci and calculate the eigenvectors and eigenvalues. Three principal components were obtained by using GCTA, and the PCA distribution map was drawn by using the R language package. The values in the axis label parentheses represent the percentage of the principal component that explains the overall variance.

Genetic Diversity and Selective Sweep Analysis
Nucleotide diversity (π) (also known as θ π, Pi), genetic diversity analysis (Tajima's D), and genetic differentiation among populations (F st ) were analyzed by the PopGenome module in the R software package. According to the phylogenetic results, P. koraiensis was divided into three groups: the Changbai Mountain group, the Great Khingan Mountains group, and the Lesser Khingan Mountains group. Because of the moderate climate and frequent gene communication with other populations in the Changbai Mountain group, we took four populations (HL, ML, LJ, and TQL) in this group as background groups. Combined with the results of F st , the MKS population of the Great Khingan Mountains group and WYL populations of the Lesser Khingan Mountains group were used as control groups to carry out selective sweep analysis. We scanned the genome in 100 kb sliding windows with a step size of 10 kb to identify significant regions based on a fixed index (F st ) and nucleotide diversity ratio (π ratio). We respectively screened the significant regions according to the top 5%, and the intersection region is the selected region. Genes in the selected region were defined as the selected gene, namely, the candidate gene. After obtaining the candidate genes in the selected region, the genes were analyzed by GO (gene ontology) (www.geneontology.org, accessed on 13 January 2023) and KEGG enrichment analysis. Then, the variations in SNPs in the coding region of the candidate genes were annotated, and the core functional mutations were speculated to provide clues for further verification.

Conclusions
In our study, phylogenetic, population structure and selective sweep analyses of different populations of P. koraiensis in northeastern China were performed by GBS-seq sequencing technology. The results of the phylogenetic analysis showed that P. koraiensis was divided into three geoclimatic zones: the Great Khingan Mountains climatic region, the Lesser Khingan Mountains climatic region, and the Changbai Mountain climatic region, with introgression in the MKS population and outgroups. Geographic isolation, climatic factors, and introgression are the driving factors for the differentiation of P. koraiensis. The MKS population at the northern edge of the distribution area and the WYL population in the mining area are the two populations with greater differentiation. The differentiation of the MKS population is consistent with the central-marginal hypothesis, and the results of selective sweep analysis show that the genes subject to selection removal in the MKS population are related mainly to flowering and photomorphogenesis, cellular response to water deficit, response to reactive oxygen species, flavonoid synthesis, and DNA repair. Climatic heterogeneity factors (e.g., sunlight, precipitation, and temperature) drive the differentiation of MKS; genes selected in the WYL population are related mainly to metal ion transport and macromolecule biosynthesis, such as flavonoids, phytohormone signaling, and DNA repair, and the WYL population may be subject to heavy metal stress to produce differentiation. Our results provide insights into adaptive evolution molecular mechanisms of genus spruce and contribute to spruce resistance improvement via genomic-assisted breeding.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/plants12061266/s1, Table S1: List of selected genes in P. koraiensis MKS (Mengkeshan) population; Table S2: List of selected genes in P. koraiensis WYL (Wuyiling) population; Table S3: GO categories of selected genes in P. koraiensis MKS (Mengkeshan) population; Table S4: GO categories of selected genes in P. koraiensis WYL (Wuyiling) population. Funding: This study was supported by the "Study of Phylogeography and Quaternary glacial refuge in Picea koraiensis (31500540)" and the "study on Phylogeny and Phylogeography in Keteleeria (32271694)" funding by the National Natural Science Foundation of China.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The sequence reads associated with the paper are stored at the NCBI Sequence Read Archive in the BioProject PRJNA664562.