Abstract The taxonomic classification of Picea meyeri and P. mongolica has long been controversial. To investigate the genetic relatedness, evolutionary history, and population history dynamics of these species, genotyping-by-sequencing (GBS) technology was utilized to acquire whole-genome single nucleotide polymorphism (SNP) markers, which were subsequently used to assess population structure, population dynamics, and adaptive differentiation. Phylogenetic and population structural analyses at the genomic level indicated that although the ancestor of P. mongolica was a hybrid of P. meyeri and P. koraiensis , P. mongolica is an independent Picea species. Additionally, P. mongolica is more closely related to P. meyeri than to P. koraiensis , which is consistent with its geographic distribution. There were up to eight instances of interspecific and intraspecific gene flow between P. meyeri and P. mongolica . The P. meyeri and P. mongolica effective population sizes generally decreased, and Maxent modeling revealed that from the Last Glacial Maximum (LGM) to the present, their habitat areas decreased initially and then increased. However, under future climate scenarios, the habitat areas of both species were projected to decrease, especially under high-emission scenarios, which would place P. mongolica at risk of extinction and in urgent need of protection. Local adaptation has promoted differentiation between P. meyeri and P. mongolica . Genotype‒environment association analysis revealed 96,543 SNPs associated with environmental factors, mainly related to plant adaptations to moisture and temperature. Selective sweeps revealed that the selected genes among P. meyeri , P. mongolica and P. koraiensis are primarily associated in vascular plants with flowering, fruit development, and stress resistance. This research enhances our understanding of Picea species classification and provides a basis for future genetic improvement and species conservation efforts.
