Abstract Wild boars inhabit diverse climates, including frigid regions like Siberia, but their migration history and cold adaptation mechanisms into high latitudes remain poorly understood. We constructed the most comprehensive wild boar whole-genome variant dataset to date, comprising 124 samples from tropical to frigid zones, among which 47 Russian, 8 South Chinese and 3 Vietnamese wild boars were newly supplemented. We also gathered 75 high-quality RNA-seq datasets from 10 tissues of 6 wild boars from Russia and 6 from southern China. Demographic analysis revealed the appearance of Russian wild boars in Far East of Asia (RUA) and Europe (RUE) after the last glacial maximum till ~ 10 thousand years ago. Recent gene flow (<100 years) from RUA to RUE reflects human-mediated introductions. Cold-region wild boars exhibit strong selection signatures indicative of genetic adaptation to cold climates. Further pathway and transcriptomic analyses reveal a novel cold resistance mechanism centered on enhanced vitamin A metabolism and catalysis, involving the reuse of UGT2B31 and rhythm regulation by ANGPTL8, RLN3 and ZBTB20. This may compensate for the pig’s lack of brown fat/UCP1 thermogenesis. These findings provide new insights into the molecular basis of cold adaptation and improve our understanding of Eurasian wild boar migration history. Subject terms: Population genetics, Agricultural genetics __________________________________________________________________ A study explores the genomic and transcriptomic mechanisms underlying vitamin A-induced thermogenesis in wild boars, revealing gene reuse as a cold adaptation strategy and providing insights into evolutionary adaptations to cold environments. Introduction Wild boar (Sus scrofa) is unique among ungulates for its practice of constructing thermoprotective nests, a behavior that provides a warm microenvironment crucial for the survival of its litter during and immediately after birth^[38]1. Unlike many mammals, piglets lack brown adipose tissue (BAT) and depend primarily on shivering for thermoregulation^[39]2. This thermoregulatory limitation is likely linked to the evolutionary history of Sus scrofa, which originated from Island South East Asia (ISEA) with a tropical or subtropical climate^[40]3–[41]6, and diverged from other pig genera approximately 1.4 million years ago (Mya)^[42]7. The disruption of the essential thermogenic gene UCP1 in the suid lineage, which encodes a brown adipocyte-specific mitochondrial uncoupling protein crucial for regulating body temperature, likely contributed to this thermoregulatory incapacity. Over long-term colonization in warm environments, the UCP1 gene underwent little or weak selection, leading to its pseudogenization in suids^[43]1. Consequently, the absence of UCP1/BAT thermogenesis mechanism results in the observed poor thermoregulation capacity in wild boars. Despite this, wild boars have demonstrated remarkable adaptability, thriving in diverse and extreme environments, including the hypoxia high-altitude Tibet Plateau^[44]8 and the frigid Russian Far East^[45]9. This rapid expansion, spanning only about 200 thousand years, underscores their adaptability. However, the transition from tropical or subtropical climates to colder environments, where heat transfer is more rapid and temperatures are lower, would have posed significant challenges. This suggests that wild boars possess alternative mechanisms for coping with cold exposure in the absence of UCP1/BAT-dependent thermogenesis. Various alternative thermogenic strategies have been proposed for pigs, such as shivering^[46]10, UCP3-mediated heat production^[47]11 and muscle-nonshivering thermogenesis involving ATPase pump activity in the sarcoplasmic reticulum^[48]12,[49]13. Additionally, leptin has been implicated in pig thermoregulation^[50]14, although the precise mechanisms remain unclear or debated^[51]13,[52]15. Previously, we performed genomic analyses to get insights on cold adaptation in domestic pigs, and identified a set of candidate cold-adaptation genes related to hair development, forebrain neuron differentiation, kidney development, energy metabolism and blood circulation^[53]11,[54]16. However, since wild boars live naturally without the interference of human-accompanied life, studying their cold adaptation mechanisms may offer deeper insights. Recent research on wild boar mitogenomes suggests that cold climate impose more stringent evolutionary pressures on the mitochondria of wild boars in colder regions; for example, Cytb gene exhibits significantly lower Ka/Ks ratios in Siberian wild boars, indicating that cytochrome b might play a critical role in cold adaptation^[55]17. Hitherto, there has been no systematic study of cold adaptation in wild boars using whole-genome analysis. Thus, a comprehensive investigation into the molecular mechanisms underlying cold adaptation in wild boars is necessary. Pigs have been extensively studied as important biomedical models for humans^[56]4. In the context of obesity and cold adaptation, pigs are more suitable models for adult humans and obese patients, as both exhibit limited levels of UCP1/BAT activity. This similarity makes pigs a more accurate model than mice (which possess abundant UCP1 and rely on BAT-driven thermogenesis) for studying human metabolism, especially in relation to obesity-induced hyperglycemia, cholesterol profiles, and energy metabolism^[57]18–[58]20. Temperature homeostasis is crucial for survival, and it is expected that thermostatic animals have evolved multiple mechanisms of thermogenesis throughout their long evolutionary history. Perspectives on pig cold adaptations are particularly valuable given their close evolutionary affinity and biological similarities to humans and may inform the development of treatments for diseases and disorders related to cold climate. Historically, high-quality genomic data on pigs were mostly limited to domesticated pigs and wild boars from China, Southeast Asia, and Europe. However, there has been a significant gap in genomic data for wild boars inhabiting the cold regions of northern Asia, particularly Siberia. Consequently, our knowledge regarding the evolutionary history and cold adaptation mechanisms of wild boars in these regions remains limited. The wild boar populations in North Asia represent a crucial piece of the puzzle in understanding the migration history of pigs. It remains unclear whether Russia, located in North Asia, served as a key corridor linking the spread of wild boars from Indonesian islands to Asia and Europe, that is, whether some Asian pigs may have dispersed to Europe via Russia. In this study, we analyzed genomic and transcriptomic data from Russian wild boars residing in extreme cold regions, including 42 wild boars from the Asian Far East part of Russia (RUA) and 5 wild boars from the European part of Russia (RUE). We also collected tissue samples from 10 different tissues of 6 Far East Russian wild boars and 6 South Chinese wild boars for transcriptome sequencing. Our objectives were to address the following questions: when did Russian wild boar populations originate? What is the population structure of RUA? Did Russia serve as an alternative route for the spread of Asian wild boars to Europe, or was there gene flow between RUA and other pig populations? What is the underlying mechanism of cold adaptation in these populations? Are their cold adaptation mechanisms consistent with those of other high-latitude wild boars? This study enhanced our understanding of the evolutionary history of pigs, deepened our understanding of mammalian cold adaptation mechanisms, and highlighted the potential of using pigs as models for cold tolerance. Results High-quality genetic data of global wild boars In this study, we collected 124 high-quality whole-genome resequencing data from global wild boars across the different climatic zones, ranging from tropical to extreme cold environments (Fig. [59]1a). Specifically, we obtained new genetic samples from 47 Russian wild boars, 8 South Chinese wild boars, and 3 Vietnam wild boars for whole-genome sequencing and analysis (Supplementary Data [60]1). To create a comprehensive dataset of wild boars, we supplemented publicly available whole genomes including 11 South Chinese wild boars (SCW), 3 North Chinese wild boars (NCW), 10 Korean wild boars (KWB), 1 Japanese wild boar (JWB), 26 European wild boars across North and South Europe, 2 Sumatran wild boars (SMW), 1 Armenian wild boar (AMW) and 12 outgroups (Supplementary Data [61]1). All 124 genomes were aligned, genotyped, and annotated using the Duroc pig reference genome (Sus scrofa 11.1). A total of 45,013,177 putative single-nucleotide polymorphisms (SNPs) were yielded and passed the filtering criteria across the 124 genomes for subsequent analyses. To assess the quality of the genetic variants, we calculated the transition-to-transversion (Ts/Tv) ratio, a commonly used metric to assess variant quality^[62]21, which was found to be 2.49 for the global wild boar SNPs in our study. This high Ts/Tv ratio suggests that the genetic variants identified in our study are of high quality. Fig. 1. Sample distribution and demographic history of Eurasian pigs. [63]Fig. 1 [64]Open in a new tab a Map showing the geographic distribution of the sampling localities of all wild boars across Eurasia and outgroup covered by our study. Terrestrial surface temperatures averaged from 1960 to 1990 A.D. are shown as a heat map with scale bars in °C. Copyright © Esri. All rights reserved. b Maximum-likelihood tree depicting the evolutionary relationships among Eurasian wild boars and outgroup. The numbers at the major branches are bootstrap values. c ADMIXTURE analysis with K = 2, 4 and 8. Colors in each column represent ancestry proportion. d Relationships among Eurasian pigs inferred using TreeMix. e Heat map of the f-branch metric for Eurasian wild boars and outgroup. The used guide tree is shown along the x- and y-axes (in “laddered” form along the y-axis). The matrix shows the inferred f-branch metric, reflecting excess allele sharing between the branch of the “laddered” tree on the y-axis (relative to its sister branch) and the branches defined on the x-axis. The asterisks denote significance with Z-score >3. f The demographic model inferred with momi2. See Supplementary Data [65]1 for full names of abbreviation. AWP, Phacochoerus africanus; SMB, Sumatran wild boar; VWB, Vietnamese wild boar; SCW, South Chinese wild boar; KWB, Korean wild boar; JPW, Japanese wild boar; NCW, North Chinese wild boar; RUA, Russian wild boar in Asian part; RUE, Russian wild boar in European part; SEW, South European wild boar; AMW, Armenian wild boar; NEW, North European wild boar; UW, Ukraine wild boar; SGEW, wild boar from Samos, Greece. Independent and rapid spread of European and Asian Wild Boars from South Asia To decipher genetic relationships among Eurasian wild boar populations, we constructed a maximum likelihood (ML) tree using the SNPs on 4-fold degenerate (4DTv) sites and conducted principal component analysis (PCA) using the whole-genome wild boar dataset (Supplementary Fig. [66]1). Our results exhibited an overall phylogenetic relationship of Eurasian wild boars and clearly positioned Russian wild boars from both European and Far East Asian regions within these relationships (Fig. [67]1b). Consistent with previous studies^[68]5, our phylogenetic tree showed that Sumatran wild boars are basal to both Asian and European wild boars. Except for one Vietnamese wild boar, European and Asian wild boars clustered separately. Within the Asian group, the phylogenic relationships suggested a northward expansion, with Vietnamese wild boars at the base of the Asian clade, followed by South Chinese wild boars (SCW), Japanese wild boars (JWB), Korean wild boars (KWB), North Chinese wild boars (NCW), and finally Russian wild boars from Far East Asia (RUA). The genetic divergence between southern and northern Asian wild boars was strongly supported (100% bootstrap). RUA clustered at the end of the Asian group, indicating that they were the most recently formed among the Asian populations. Furthermore, RUA could be divided into two smaller genetic clusters: one comprising individuals from Terney and Arharinsky districts (named as RUA1) and another from Blagoveshensk district (named as RUA2), with RUA1 showing close genetic affinity to NCW (Fig. [69]1b). In contrast, the phylogenetic relationships among European wild boars did not suggest a straightforward northward invasion route. Notably, Ukrainian wild boar (UW) and RUE were located on deeper split branches than the other European populations (Fig. [70]1b). This positioning could be explained by two hypotheses: (1) among the earliest formed European wild boar populations, or (2) they were the result of admixture between European wild boars and Asian wild boar ancestry, skewing their phylogenetic positions towards the Asian group. To distinguish between these two hypotheses, we conducted admixture and migration analyses. We performed individual ancestry inference to assess the population structure using ADMIXTURE^[71]22 for population structure inference. The optimal K value revealed by ADMIXTURE was four (Fig. [72]1c, Supplementary Fig. [73]2), at which all individuals grouped into four genetic clusters corresponding to their geographic distributions (outgroups in ISEA, wild boars in South Asia, ones in North Asia and ones in European; Fig. [74]1c). Notably, the ancestral component of Asian wild boars was only detected in UW and RUE located near the basal of European wild boars. UW contained ancestry of South Asian wild boars, while RUE’s Asian ancestry originated from North Asian wild boars. As K increased, we identified RUE’s Asian ancestry as originating from RUA1 in Terney and Arharinsky districts (Fig. [75]1c, Supplementary Fig. [76]2). TreeMix results with optimal migration events (equal to 6 times) also supported RUE’s Asian ancestry was from RUA1, and that the Asian ancestry of UW was likely derived from the ancestor of South Asian wild boars (Fig. [77]1d, Supplementary Fig. [78]3). Additionally, there was gene flow from SCW to VWB. All of these migration events were corroborated by the F-branch analyses^[79]23,[80]24, which confirmed three gene flows between RUA1 and RUE, between UW and Asian wild boars, and between SCW and VWB (Fig. [81]1e). Notably, the strongest gene flow (Fig. [82]1e; dark red) was detected between RUE and ancestor of European wild boars in the F-branch analyses. To further investigate the demographic histories of these wild boar populations, we inferred long-term effective population size (N[e]) dynamics for each population and using pairwise sequential Markovian coalescent (PSMC)^[83]25 and MSMC2^[84]26 methods. The N[e] trends of RUA and RUE (Supplementary Fig. [85]4a) were consistent with those previously reported for Asian and European wild boars^[86]7. MSMC splitting time exhibited that divergence between RUA and NCW occurred ~8876 years ago while divergence time between RUE and European wild boars was dated to ~ 8122 years ago (Supplementary Fig. [87]4b). Both split times were the most recent among all population pairs. However, past studies have indicated that significant introgression can bias divergence time estimates when using relative cross coalescence rates, especially when mixture proportions exceed 30%^[88]27,[89]28. To address this, we combined phylogeny, genetic admixture, and N[e] estimation and modeled the demographic history of Eurasian wild boars using site-frequency spectrum (SFS) via momi2^[90]29. We varied different models (such as: without migration, times pulse of migration, and direction of migration) to seek the best one with lowest AIC value (see Materials and Methods and Supplementary Fig. [91]5). The optimal model (Fig. [92]1f) suggested that ancestors of RUA and NCW diverged ~ 9855 years ago (95% CI, 8.4 to 12.0 Kya), while RUE split from EUW ~ 7041 years ago (95% CI, 6.2–8.9 Kya). Assuming pulse-like migration, two significant gene flow events were inferred from European wild boars to RUE occurring 3548 and 45 years ago with migration rates of ~40.3% and 29.0%, respectively. Additionally, RUE received ~3.4% gene flow from RUA 45 years ago. Other notable gene flow events included ~40.0% from SCW to VWB 16,906 years ago and ~10.9% from European wild boars to SCW 7562 years ago. These results confirm that admixture influenced the phylogenetic positioning of UW and RUE, and support the hypothesis that European wild boars followed a south-to-north invasion route. In summary, our results provided a comprehensive overview of the evolutionary history of wild boars to date, demonstrating that both European and Asian groups independently and rapidly radiated from South Asia to Europe and other parts of Asia within the past 200 thousand years. The evolutionary position of Russian wild boars, situated at the periphery of the wild boar phylogeny, is consistent with previous findings that wild boars reached the Eurasian continent around 219 Kya. Subsequently, some migrated to Asia while others moved towards Europe, eventually spreading northward and adapting to various climates, including the severe cold regions of present-day Russia. RUA and RUE appeared around 10 Kya and 7 Kya in Asian and European parts of Russia, respectively. This time was intersection with that of human hunting and pig domestication, as well as an increase in temperature after the Last Glacial Maximum. These may have been main factors driving wild boars’ northward migration. Additionally, substantial gene flow from European wild boars to RUE (~29.0%) and from RUA to RUE (~3.4%) occurred within the past 100 years. Reuse of UGT2B31 in Asian and European lineages We applied a composite selection signal (CSS) test^[93]30, a weighted index of selection signatures from multiple estimates (π-ratios, F[ST], and XP-CLR) to identify candidate genomic regions (50 kb in length) with robust signatures of selection in RUA population, where South Asian wild boars (including SCW and VWB) was used as a control population (Fig. [94]2; see Materials and Methods). The top 1% regions of CSS (P < 0.035) were regarded as candidate regions across the genome that had significant evidence of positive selection^[95]31. This approach allowed for identification of a total of 340 candidate genes harbored in 452 loci (Supplementary Data [96]2). At the pathway level, no significant GO term was enriched by these positively selective genes. However, several significant metabolism-related KEGG pathways, such as retinol metabolism and metabolism of xenobiotics by cytochrome P450, were associated with the candidate genes bearing selective signals (Table [97]1). These results suggested that natural selection might have acted on wild boars to enhance or modify their metabolic processes in response to cold climates. Additionally, the thyroid hormone signaling pathway, which plays a crucial role in energy homeostasis and cold adaptation^[98]32, was significantly enriched (P = 0.03). Fig. 2. Selective sweep regions in RUA. [99]Fig. 2 [100]Open in a new tab a Genomic landscape of selection signal in RUA. b Signals of positive selection on UGT2B31. UGT2B31-located region showed a lower diversity (π and π-ratio) and higher differentiation (iSAFE) compared with SCW. Orange dots depict the SNPs with highest iSAFE score around UGT2B31 in RUA, and red dot indicates a missense mutation. c The multi-species alignment of the nonsynonymous substitution encoded in the UGT2B31 gene (red dot in b). The amino acid sequences in northern wild boars around the substitution shown here are same as those in southern wild boars, except amino acid at site of 314. Dots indicate identity with the master sequence. d The geographical distribution of the ancestral (red) and derived (blue) alleles at chr8:66,317,760 in wild boars. A straight-line fitting indicated a significant negative correlation between the derived allele frequency and temperature in wild boars. e Homology modeling of protein structure of UGT2B31 containing ancestral and derived allele, respectively. f The geographical distribution of the ancestral (red) and derived (blue) alleles at chr8:66,317,760 in domestic pigs. A straight-line fitting indicated no significant correlation between the derived allele frequency and temperature in domestic pigs. EULD, European local domestic pig; PT, Pietrain; LD, Landrace; LW, Large White; DU, Duroc; WDU, White Duroc; MIN, Min pig; HT, Hetao pig; BAM, Bamei pig; LWU, Laiwu pig; SCT, Sichuan Tibetan pig; YNT, Yunnan Tibetan pig; BS, Baoshan pig; NJ, Neijiang pig; BMX, Bamaxiang pig; WZS, Wuzhishan pig; EHL, Erhualian pig; JH, Jinhua pig. Terrestrial surface temperatures averaged from 1960 to 1990 A.D. are shown as a heat map with scale bars in °C. Copyright © Esri. All rights reserved. Table 1. KEGG pathway enrichment analysis with positively selected genes detected by CSS ID Term Adjusted p value Gene Count Gene name ssc00830 Retinol metabolism 4.48 × 10^–8 13 UGT2B31, LOC100626199, ALDH1A3, CYP3A22, HSD17B6, LOC106510546, LOC100623504, LOC100624788, LOC100512656, RDH16, LOC100516628, LOC100624700, LOC100515394 ssc00040 Pentose and glucuronate interconversions 1.61 × 10^–6 8 UGT2B31, SORD, CRPPA, LOC100623504, LOC100624788, LOC100516628, LOC100624700, LOC100515394 ssc00983 Drug metabolism - other enzymes 2.13 × 10^–6 11 UGT2B31, GSTA4, UPB1, LOC100510917, DUT, LOC100623504, LOC100624788, LOC100516628, LOC100624700, LOC100526118, LOC100515394 ssc00980 Metabolism of xenobiotics by cytochrome P450 4.04 × 10^–6 10 UGT2B31, GSTA4, LOC100510917, LOC106510546, LOC100623504, LOC100624788, LOC100516628, LOC100624700, LOC100526118, LOC100515394 ssc05204 Chemical carcinogenesis - DNA adducts 4.04 × 10^–6 10 UGT2B31, GSTA4, CYP3A22, LOC100510917, LOC100623504, LOC100624788, LOC100516628, LOC100624700, LOC100526118, LOC100515394 ssc00140 Steroid hormone biosynthesis 6.51 × 10^–6 10 UGT2B31, SULT2B1, CYP3A22, LOC100620470, LOC100623504, LOC100624788, SULT1E1, LOC100516628, LOC100624700, LOC100515394 ssc04976 Bile secretion 1.32 × 10^–5 11 UGT2B31, AQP9, ATP1B3, ABCB4, SLC22A7, LOC100522267, LOC100623504, LOC100624788, LOC100516628, LOC100624700, LOC100515394 ssc00982 Drug metabolism - cytochrome P450 1.32 × 10^–5 9 UGT2B31, GSTA4, LOC100510917, LOC100623504, LOC100624788, LOC100516628, LOC100624700, LOC100526118, LOC100515394 ssc05207 Chemical carcinogenesis - receptor activation 2.31 × 10^–5 16 UGT2B31, PAQR8, GSTA4, MTOR, CREB3L3, CYP3A22, LOC100510917, MAP2K2, EPHX3, PRKCB, LOC100623504, LOC100624788, LOC100516628, LOC100624700, LOC100526118, LOC100515394 ssc00053 Ascorbate and aldarate metabolism 8.97 × 10^–5 6 UGT2B31, LOC100623504, LOC100624788, LOC100516628, LOC100624700, LOC100515394 ssc00860 Porphyrin metabolism 9.83 × 10^–4 6 UGT2B31, LOC100623504, LOC100624788, LOC100516628, LOC100624700, LOC100515394 ssc04919 Thyroid hormone signaling pathway 3.35 × 10^–2 8 SLCO1C1, MTOR, PLCE1, ATP1B3, ITGAV, NOTCH3, MAP2K2, PRKCB [101]Open in a new tab Notably, the gene UGT2B31 was involved in nearly all significant enriched KEGG pathways, except for the thyroid hormone signaling pathway. UGT2B31 is known to catalyze the binding of lipophilic substrates with glucuronic acid, thereby increasing the water solubility of metabolites. For example, it converts retinoic acid (RA) to retinoyl-β-glucuronide (RAG) in the retinol metabolism pathway proceeding in liver^[102]33. We detected a strong selective sweep spanning a 1.25-Mb region (65.75–67.00 Mb) harboring UGT2B31 (Fig. [103]2b). Further analysis using iSAFE^[104]34 identified the SNP hotspots favored by ongoing selective sweep within this region. The results showed that the SNPs with highest score were distributed as a horizontal line across UGT2B31, suggesting strong linkage disequilibrium in the region. Notably, a missense mutation within the UGT2B31 gene (chr8:66,317,760 A/C: UGT2B31-Ser314Arg) was observed in RUA (Fig. [105]2b, red dot; Fig. [106]2c, Supplementary Fig. [107]6). Interestingly, the genotype at this locus showed a significant correlation with temperature (Fig. [108]2d). Detailly, as the temperature decreased, the derived allele became more frequent in Eurasian wild boars, tending to be fixed when temperatures dropped below 15 °C (Fig. [109]2d). This pattern was observed in both Asian and European wild boar lineages. The simultaneous fixation of the derived allele in these independent lineages cannot be explained by the recent migration between Asian and European wild boars, because it required extensive gene flow between them, but such a large-scale introgression event has not been suggested in the population history analysis (Fig. [110]1). Instead, the selection pressure from cold environments likely drove the simultaneously gradual fixation of the derived allele in both lineages. This is a phenomenon of gene reuse, where the same underlying genes are under repeated selection among independently radiating populations. It emphasizes the key role of UGT2B31 in cold adaptation. We further performed protein structure predictions at [111]https://robetta.bakerlab.org to contrast two protein structure containing ancestral and derived genotype. Although the overall structure of the protein did not exhibit significant changes, the substitution of Ser with Arg at the target amino acid, located in the “hook” structure (Fig. [112]2e), may have enhanced substrate capture and binding capacity, potentially increasing the efficiency of the UGT2B31-encoded protein. To extend our findings, we included domestic pig breeds from our previous study^[113]35. We examined the distribution of the missense mutation (chr8:66,317,760) in Eurasian domestic pigs. The results showed that the derived allele frequency in domestic pigs followed a similar trend with temperature as in wild boars, though it did not reach statistical significance (Fig. [114]2f). This may be due to the influence of human activities, such as the rearing of livestock in indoor environments and human-mediated gene flow, which could have disrupted the natural selection signals at this locus. As a result, when studying the genetic mechanism of cold adaptation in domestic pigs, we could not dig out the positive signal of this locus in our previous study^[115]16. Additional selected genes associated with thermoregulation In our analysis, we identified several genes located in top selective sweep regions that are associated with thermoregulation and other related physiological processes. Notable among these are ANGPTL8, RLN3, and LOC100516957, which exhibited composite selection signal (CSS) scores of 1.66, 1.91, and 2.28, respectively (Fig. [116]2a, Supplementary Data [117]2). According to their annotations in the NCBI Gene database, these genes are functionally relevant to circadian rhythms and metabolic regulation. ANGPTL8 (also known as betatrophin and Lipasin) belongs to the angiopoietin-like protein family and is a circulating protein highly enriched in liver and adipose tissue^[118]36,[119]37. ANGPTL8 plays a role in mediating the circadian rhythms of liver, the body’s largest metabolic organ, which responds sensitively to food signals and secretes hepatokines, thereby regulating metabolic and circadian processes^[120]37. Interestingly, ANGPTL8 can also be thermoregulated; in mice, cold exposure induces its expression in BAT more than threefold, suggesting its distinct physiological roles in mediating thermogenesis and energy homeostasis^[121]38. In addition, ANGPTL8 is involved in lipid metabolism and exhibits a negative relationship with plasma high density lipoprotein-cholesterol^[122]39. RLN3 is another gene related to the regulation of circadian rhythms^[123]40 and is known to simulate food intake^[124]41. Similarly, LOC100516957 is highly expressed in neurons that govern food intake, further linking it to energy balance and metabolic processes^[125]42. There are also some positively selected gene associated with mitochondrial function. For example, TIMM21, located on chromosome one with CSS score of 1.88, is involved in the assembly of mitochondrial cytochrome c oxidase and encodes components of mitochondrial membrane^[126]43. MTX2 (CSS score: 1.72) encodes Metaxin-2, an outer mitochondrial membrane protein involved in the transport of proteins into mitochondrion^[127]44. Additionally, NDUFAF2 (CSS score: 1.51) and NDUFB1 (CSS score: 1.89) are crucial for the assembly of mitochondrial respiratory chain complex I^[128]45. Another gene, ACSS1 (CSS score: 2.14), encodes mitochondrial acetyl-CoA synthetase, which plays a specific and unique role in thermogenesis during fasting in mice^[129]46. Furthermore, the gene LDHB, located in a region under positive selection with a CSS score of 2.14, encodes a key glycolytic enzyme responsible for converting lactate and NAD^+ to pyruvate, NADH and H^+^[130]47. Pyruvate produced is then utilized in the tricarboxylic acid (TCA) cycle, leading to increased heat production^[131]48. Another gene under strong positive selection is IGF1R, with a CSS score of 2.06. This gene encodes the insulin-like growth factor 1 receptor, a receptor-tyrosine kinase that plays a critical role in signaling pathways important for cell survival and proliferation^[132]49. IGF1R has been implicated in the regulation of body across various species, including dogs^[133]50, cattle^[134]49, mouse^[135]51 and human^[136]52. Genes and pathways responding to cold at the transcriptomic level To further elucidate the mechanisms by which wild boars respond to severe cold climates, we generated and analyzed 75 transcriptomic datasets derived from 10 organs or tissues, including the hypothalamus, liver, heart, muscle, skin, kidney, spleen, stomach, bone and subcutaneous fat, across 12 wild boars (6 Russian wild boars and 6 South Chinese wild boars) (Fig. [137]3, Supplementary Data [138]3). Hierarchical cluster analyses revealed distinct transcriptomic profiles between the RUA and SCW samples in nearly all organs or tissues, with the exception of the stomach (Supplementary Fig. [139]7). A total of 22,063 Sus scrofa reference genes downloaded from the Ensembl database (version 108). Among these genes, 3610 genes (1819 upregulated, 1791 downregulated), 1858 genes (1060 upregulated, 798 downregulated), 2302 genes (1232 upregulated, 1070 downregulated), 4619 genes (2223 upregulated, 2396 downregulated), 3576 genes (1841 upregulated, 1735 downregulated), 2348 genes (1404 upregulated, 944 downregulated), 1290 genes (560 upregulated, 730 downregulated), 2738 genes (1694 upregulated, 1044 downregulated), 6893 genes (3727 upregulated, 3166 downregulated) and 3382 genes (1991 upregulated, 1391 downregulated) were detected as being differentially expressed between RUA and SCW in hypothalamus, liver, heart, muscle, skin, kidney, spleen, stomach, bone and subcutaneous fat, respectively (Supplementary Data [140]4–[141]13). Fig. 3. Results of RNA-seq related analyses. [142]Fig. 3 [143]Open in a new tab a Results of enrichment analysis with upregulated genes of 10 tissues in RUA. b Upset plot of tissue-specific and tissue-shared unique upregulated genes. c ClueGO network for the unique upregulated genes in liver. d Tissues with significantly upregulated expression of COX2. e Tissues with significantly upregulated expression of ZBTB20. (*: P < 0.05, **: P < 0.01, ***: P < 0.01). KEGG and GO enrichment analysis was conducted for the up-regulated genes in all RUA tissues. The results revealed that the retinol metabolism pathway (KEGG: ssc00830) was significantly enriched in seven organs or tissues, including the hypothalamus, liver, heart, skin, stomach, bone, and subcutaneous fat (Fig. [144]3a; see Supplementary Data [145]14–[146]23 for detailed information). The thermogenesis pathway (KEGG: ssc04714) was enriched in four organs, including heart, muscle, spleen and stomach. Additionally, pathways such as oxidative phosphorylation (KEGG: ssc00190), steroid hormone biosynthesis (KEGG: ssc00140) and ribosome (KEGG: 03010) were significantly enriched in multiple organs or tissues in RUA (Fig. [147]3a). Beyond these shared pathways, several organ-specific pathways were observed. For example, the hippo signaling pathway was enriched in bone tissue, while the citrate cycle (TCA) and pyruvate metabolism pathways were enriched in muscle. In the hypothalamus, pathways such as circadian entrainment, calcium signaling, and thyroid hormone synthesis were significantly enriched. These pathways are involved in rhythms, metabolism, and growth, and are likely utilized by RUA to respond to cold stress. Subsequently, we investigated the upregulated genes that were either unique to each organ or shared across different tissue combinations (Supplementary Data [148]24). The retinol metabolism remained enriched by genes (UGT2B31, LRAT and CYP26A1) unique to liver (Fig. [149]3b). Network analyses via GlueGO suggested that retinol metabolism could indirectly regulate the mitochondrial membrane potential (GO:0051881) (Fig. [150]3c). Additionally, the gene ANGPTL8, which regulates food intake and rhythm, was uniquely and significantly upregulated in the liver. UGDH, which catalyzes the oxidation of UDP-glucose to UDP-glucuronate (used by UGT2B enzymes in glucuronidation)^[151]53, was also upregulated in the liver. Genes involved in RA or retinol transport, such as TTR and ALB, were significantly upregulated in eight organs or tissues in RUA (Fig. [152]3b). Notably, the COX2 gene was abundantly expressed in visceral organs and subcutaneous fat of RUA, with expression fragments per kilobase of exon per million fragments mapped (FPKM) values ranging from 10,000 to 90,000, significantly higher than in SCW (Fig. [153]3d). This suggests that highly active systemic mitochondrial synthesis and mobilization are critical for coping with cold temperatures. Intriguingly, ZBTB20 was significantly upregulated in all 10 organs or tissues (Fig. [154]3e). This gene is closely related to rhythm regulation, and its loss has been shown to impair circadian output and leads to unimodal behavioral rhythms^[155]54. This finding strongly indicates the key role of rhythm regulation in maintaining body temperature balance in cold environments. Wild boars from Far East are generally larger than most South and South-East Asian pigs and exhibit less variability in size^[156]55. We noted that IGF1R was down-regulated in almost all RUA tissues compared to SCW (Supplementary Fig. [157]8), consistent with the expression pattern in other mammals where reduced expression of IGF1R contributes to larger body size within same species^[158]49–[159]52. A larger body size is advantageous in colder climates, as predicted by Bergmann’s Rule, because it allows for a lower mass-specific rate of energy expenditure and enhances an organism’s ability to conserve heat^[160]56. Discussion In this study, we systematically explored the origin and demographic history of Russian wild boars in the broader context of the evolution of Eurasian wild boars. A recent study speculated that Russia might have served as an alternative route for wild boars to spread from Asia to Europe due to its vast geographical spanning^[161]57. However, our findings suggest that this hypothesis is not supported by the divergence times of Russian wild boars. In our previous study, we found that the formation of North China wild boars coincided with the timing of human hunting and domestication activities, and occurred much more recently than the emergence of European wild boars^[162]7. This indicates that large-scale human activities may have driven wild boars to migrate from warmer environments to colder ones. In this study, we inferred the formation time of RUA was about 10 Kya, less than that of European wild boars, which further consolidated that Russia could not have been the route of wild boar expansion. However, we detected the genetic ancestry from RUA in RUE, which at first seemed contradictory. After detailed historical reconstruction of Eurasian wild boars, we found that the genetic components of RUA present in RUE were the result of recent introgression events over the past 100 years, more likely mediated by human activities. Historical records showed that Russian wild boars in Europe were nearly driven to extinction by large-scale hunting during the 15th century^[163]58. In the 1950s, the former Soviet Union initiated a population recovery program, which involved introducing wild boars from peripheral regions, including the Far East, into European parts of Russia^[164]59. This event not only coincided with the timing of gene flow from RUA to RUE (Fig. [165]1f), but also aligned with the results suggested by F-branch analysis that a variety of European wild boars were mixed into the RUE and the recent European wild boar migration into the RUE detected through demographic analysis. So far, by applying high-quality genomic data from Russian wild boars, we have completed a crucial piece of puzzle regarding the spread of wild boars in Asia after their departure from ISEA. Our findings confirmed that Russia was unlikely to have been the route for wild boars spread from Asia to Europe, and that the ancestral components of RUA found in RUE were the result of recent human-mediated introductions. Selection signals and differentially expressed genes suggest that retinol (vitamin A) metabolism pathway plays a critical role in the cold adaptation of wild boars. Vitamin A, a fat-soluble vitamins, generally exists in three main forms (retinol, retinal and retinoic acid), collectively known as retinoids^[166]60. Hereafter, we refer to retinoids synonymously as vitamin A. RA, a biologically active form of vitamin A (retinoids), has all-trans-RA (ATRA) as its main derivative. Literature-mining research indicates that vitamin A is an important winter nutrient. An increased requirement for vitamin A in cold conditions was first indicated by an experiment showing that vitamin A-deficient rats had a reduced survival time when exposed to cold. This study found that at least 20 times more retinoic acid was required to maintain growth and survival in cold conditions compared to 25 °C^[167]61. Subsequent studies gradually revealed three pivotal thermogenesis-related function of RA: overall energy balance, fatty acid oxidation and rhythms^[168]62–[169]65. Specifically, experiments demonstrated that RA treatment increased expression of genes related to oxidative metabolism, fatty acid oxidation, and thermogenesis in white adipose tissue (WAT) depots^[170]64,[171]66. Additionally, RA activates BAT and enhances lipid oxidation capacity in tissues such as skeletal muscle^[172]67 and the liver^[173]68. RA also affects mitochondrial biogenesis and function. Previous transcriptomic analysis indicated that RA impacted mitochondria in adipocytes, leading to increased oxidative phosphorylation capacity, oxygen consumption rate and mitochondrial number, particularly in WAT^[174]69–[175]71. Moreover, RA plays a homeostatic role as a regulator of biological rhythms in the central nervous system^[176]72, influencing the body’s response to both the day-night cycle and seasonal physiological regulation. A key function of rhythms is to anticipate favorable times of year for seasonal behaviors such as reproduction and hibernation, and nutrients are known to inform the clock about the quality and availability of food in the environment^[177]72. Rhythm- related genes, e.g. ANGPTL8, RLN3, LOC100516957 and ZBTB20, detected in this study align well with the previous research. Notably, systemic retinoid levels are up-regulated by cold exposure. Fenzl et al. ^[178]71 detected that cold-induced increases in circulating retinol occur alongside a decreased respiratory quotient and higher lipid oxidation in humans. They further suggested that longer cold challenges could have even more pronounced effects on retinoid metabolism and biological outcomes. Conversely, retinoid deficiency has been reported to result in decreased expression of mitochondrial genes in rats^[179]73. In our study, we found that the mitochondrial genes were highly expressed in RUA. For example, COX2, involved in the mitochondrial electron transport chain driving oxidative phosphorylation, was differentially highly expressed in seven RUA tissues (Fig. [180]3b), suggesting an increase in retinoids under cold conditions. In addition, we also performed a mitochondrial quantification on wild boar. We further found that the mitochondrial content of RUA was significantly higher than that of SCW (Supplementary Fig. [181]9). The high expression of COX1/2/3 genes, which are only translated on mitochondrial ribosomes^[182]74, indicates a substantial requirement for ribosomes to translate these genes, which explains observed enrichment of ribosome pathway in all the experimental tissues. Thus, we propose a plausible hypothesis that cold exposure leads to an increase of retinoids, which, in turn, enhances the number and efficiency of mitochondria, increases higher lipid utilization, boosts cellular respiration and thermogenic capacity across multiple organs. While cold exposure can induce an increase in retinoids, excessive accumulated vitamin A can be toxic for the body^[183]75. In present study, we observed a physiological coping strategy of wild boars that helps maintain high levels of vitamin A. Our results regarding selection and expression demonstrated that, after mediation by a series of genes such as ALDH1A1/2/3, AOX1 and HSD17B6, retinol was ultimately converted into RAG by UGT2B31 (Fig. [184]4), which was simultaneously under positive selection and overexpressed in the livers of Russian wild boars (Figs. [185]2, [186]3). These results also suggested that RAG might serve a physiologically significant role in vitamin A function. Although ATRA is the active form of vitamin A^[187]76, numerous studies have indicated that the biological activity of RAG is similar to that of ATRA, but without the toxic side effects of RA^[188]77,[189]78. Interestingly, even a 10-fold excess of RAG was not toxic to the cells compared to RA. It has been shown that RAG, like RA, can support the growth of vitamin A-deficient rats similarly to RA, reflecting a physiologically significant role for RAG in vitamin A function^[190]79,[191]80. Furthermore, RAG is even more functionally active than RA^[192]81. In general, we found that the UGT2B31 gene in cold-adapted wild boars was under strong positive selection, and concurrently, the mutant UGT2B31 were highly expressed (Supplementary Fig. [193]10). This gene efficiently converted RA into non-toxic but functionally similar and even more active RAG, thereby increasing tolerance to retinol and enhancing oxidation efficiency and thermogenic capacity in RUA. Fig. 4. Proposed Vitamin A-dependent cold tolerance model for wild boars in cold regions. [194]Fig. 4 [195]Open in a new tab Red genes are under positive selection in RUA, blue genes exhibit upregulated expression in RUA. NDs, genes encoding NADH: ubiquinone oxidoreductase core subunit in mitochondria; COXs, COX1, COX2, COX3; ATPase, genes encoding mitochondrial membrane ATP synthase. The plot is drawn by Figdraw ([196]www.figdraw.com). Meanwhile, consistent with its important function, RAG is widely detected in vivo and found in blood and various tissues^[197]82–[198]84. After an intraperitoneal dose in rats, RAG could be detected in the blood, liver, intestine and kidney, predominating in nearly all tissues during the following 24 h, supporting a potentially important physiological role for RAG in overall vitamin A function^[199]79. How are retinoids transferred throughout the body to perform essential functions? Research suggests retinoids are transported in plasma via several mechanisms depending on their form. RA, although lacking a specific plasma protein, can bind nonspecifically to serum albumin (ALB)^[200]85, which was observed overexpressed in RUA in this study (Fig. [201]3). As a water-soluble form of RA, RAG is transported from liver and intestine, where it is mainly formed, to peripheral tissues. Alternatively, Transthyretin (TTR) is a plasma protein secreted mainly by the liver and choroid plexus^[202]86. One of pivotal role of TTR is to transport retinoids from their main storage site in the liver to target tissues. Cold stimulation in mice and humans further leads to increased circulating retinoids^[203]71. All the above studies explain why the retinoid transport genes, e.g., TTR and ALB, were significantly overexpressed in almost all sampled organs of RUA (Fig. [204]3). Notably, expression of two retinoic acid binding proteins, CRABP1/2, was very low in this study, indicating that vitamin A is not mainly transported as retinoic acid but rather as RAG, as RAG does not bind to cellular retinoic acid binding proteins (CRABPs)^[205]85. The ancestors of wild boars originated in tropical regions, and lost the important NST-related gene UCP1 approximately 20 Mya^[206]1. As they migrated to the Eurasian continent and entered North Asia ~25 Kya, they adapted the cold climate^[207]7 and formed a distinct phylogenetic cluster form their southern counterparts. In our previous study^[208]16, we explored the cold resistance of pigs using domestic breeds but found that this approach lacked the power to reveal the selection signal of UGT2B31 (Fig. [209]2f), possibly due to environmental interference from human activity. However, in wild boars, we identified a mechanism of cold adaptation involving vitamin A through the detection of the selection signals in genes related to retinol metabolism. This was particularly evident through the simultaneous reuse of UGT2B31 gene across Eurasian lineages and was further supported by expression evidence from transcriptome data (Fig. [210]3). Vitamin A is an essential winter nutrient, but excessive amounts of vitamin A and its active metabolite RA can be toxic to both the body and mitochondria^[211]87. Wild boars have evolved to use UGT2B31 to convert vitamin A into RAG, which is not only non-toxic but also more active than RA. This conversion enhances oxidation capacity, heat production capacity and fat differentiation in multiple organs and tissues, increases mitochondrial number and efficiency, and particularly improves oxidation and thermogenesis in white adipose tissue (Fig. [212]4). This adaptation is critical for pigs, which lack brown fat and rely primarily on white fat. The high expression of genes related to retinoids transport and the sharp increase in mitochondrial content indicate that wild boars increase retinoids levels in response to cold temperatures to maintain homeostasis. As a result, it is estimated that wild boars in European part of Russia lose 30–50% of body weight by the end of winter^[213]88. Moreover, thermogenesis is physiologically regulated by circadian rhythm^[214]89. During the long winter months when food is scarce, wild boars must accumulate sufficient fat reserves before winter arrives. Compared with their counterparts in warmer South Asia, Russian wild boars must also enhance their dietary intake, biological, and seasonal rhythms to cope with the challenges of survival in harsh conditions. Circadian rhythms integrate thermal, hormonal, and nutritional information through hypothalamic circuits involved in thermoregulation^[215]90,[216]91. Retinoids can influence these rhythms by mediating neuronal plasticity and neurogenesis in both the hippocampus and hypothalamus^[217]72. At the gene level, we detected a series of rhythm-regulating genes among positively selected and differentially expressed genes, such as ANGPTL8, RLN3 ZBTB20 et al. Notably, ZBTB20 was significantly expressed in the 10 tissues of RUA we collected. Notably, the enrichment of the thyroid hormone signaling pathway observed in our RNA-Seq dataset (Fig. [218]3a) may be influenced by the winter sampling conditions. Thyroid hormone signaling is crucial for metabolic regulation and thermogenesis, often upregulated in response to cold exposure. Our sampling in late December coincided with mid-winter, with temperatures ranging from –27 °C to –10 °C in Primorsky Krai, Russia, and 5–12 °C in Jiangxi, China. These environmental factors, particularly cold stress, may have impacted gene expression patterns. While our analysis focuses on genetic mechanisms, it is important to acknowledge potential environmental contributions. Future studies incorporating multi-seasonal sampling could further disentangle the effects of environmental and genetic factors on gene expression. In summary, our study provides the first evidence that vitamin A plays an essential role in cold adaptation of wild boars lacking UCP1/BAT thermogenesis, even though earlier studies have shown that vitamin A is an important nutrient for cold resistance. By reusing UGT2B31 to enhance vitamin A metabolism and regulating rhythm via ANGPTL8, RLN3 and ZBTB20 in cold-region wild boars (Fig. [219]4), we have identified a pivotal alternative mechanism of heat production that may compensate for defects in brown fat/UCP1 thermogenesis and enable wild boars to survive in cold regions. Together, our genomic, phylogenetic, and demographic analyses have uncovered the origin and genetic structure of Russian wild pigs, including the RUA and RUE populations. From a broader perspective, we have refined our understanding of pig evolution and completed the puzzle of wild boar migration history by incorporating Russian wild boar samples. The thermogenesis mechanism induced by vitamin A- including mitochondrial number, oxidative efficiency, white fat use, and rhythm regulation- represents an important alternative to BAT/UCP1 thermogenesis in wild boars that warrants further investigation. Moreover, this study suggests RAG could have potential applications in large mammals in large mammals and in humans for coping with cold climates and may play a critical role in treating related diseases. Our results provide an unprecedented genomic resource base for ongoing molecular breeding and functional research in both medicine and in agriculture. Materials and methods Sample collection and genome sequencing We collected 47 wild boar samples for whole-genome sequencing, including 42 wild boars in Far East part of Russia, 5 ones in European part of Russia, 3 Vietnam wild boars and 10 southern Chinese wild boars. The age of all wild boars ranged approximately from 1 to 3 years, and detailed records of their sex are provided in Supplementary Data [220]1. Total genome DNA was extracted and purified from blood or muscle using phenol-chloroform method. For each sample, 1–15 µg of genomic DNA was fragmented into 200–800 bp pieces using the Covaris system (Life Technologies). The resulting fragments were then processed following the Illumina DNA sample preparation protocol. This involved end-repairing, adding an A-tail, ligating paired-end adaptors, and amplifying the library via PCR with 500-bp insert sizes. Sequencing was carried out on the HiSeq 2000 platform (Illumina), generating 100-bp paired-end reads, in accordance with the manufacturer’s standard protocols. 64 previously published genomes for wild boars and outgroups were integrated into this study, resulting in a 124-sample high-quality dataset. We firstly used fastp^[221]92 to perform quality control of reads files. Then paired-end resequencing reads were mapped to the Sus scrofa reference genome (11.1)^[222]93 with BWA “BWA-MEM” algorithm^[223]94 using the default parameters. The mapped reads were subsequently processed by sorting, duplicate marking, indel realignment, and base quality recalibrating using Picard ([224]http://picard.sourceforge.net, last accessed January 5, 2021) and GATK^[225]95. A two-round procedure of SNP calling was performed using Platypus^[226]96. We applied the criteria of genotype filtration as detailly described previously when running Platypus. All SNPs were then filtered with MAF > 0.01 and SNP call rates >80%. The high-quality variants were kept for subsequent analysis. Ts/Tv ratio was calculated via VCFtools^[227]97 to validate the quality of our SNPs set. Sample distribution maps were generated using ArcGIS® software developed by the Environmental Systems Research Institute, Inc. (Esri). The base map layer, depicting terrestrial temperature averages from 1960 to 1990 A.D., was sourced from ArcAtlas™ ([228]www.arcgis.com). Phylogeny, PCA, structure and demographic analysis Maximum-likelihood tree was built using iqtree program (version: 1.6.12)^[229]98 in ASC model based on SNPs on fourfold degenerate (4Dtv) sites with optimal substitution model GTR and with 1000 bootstrap replicates. To visualize the genetic relationships among samples, we performed a PCA using package “SNPRelate” in R^[230]99 based on whole-genome SNPs dataset. ADMIXTURE v1.3.0^[231]22 was used to quantify the genomewide admixtures among all samples. Admixture was run for each possible genetic clusters (K = 2 to 9) with 10 replicates for each K. We ran CLUMPAK^[232]100 on ADMIXTURE outputs to generate consensus solutions for each K and finally identified the optimal K cluster supported by the data based on cross validation (CV) errors. Population splitting and admixture analyses were further carried out using TreeMix program^[233]101. This method requires unlinked markers, so we used SNPrelate to prune out markers with high linkage (r^2 ≥ 0.2 within 100 Kbp). We selected an appropriate number of migration edges based on the decay in variation explanation improvements with successive numbers. BITE^[234]102 was used to do 100 bootstrap replicates, and to visualize the consensus trees with bootstrap values and migration edges when doing TreeMix analyses. Another way of looking at gene flow between populations is through the D-statistics related statitics. Herein, we implemented the f-branch test^[235]23,[236]24 to uncover correlations between results that may indicate ancestral gene flow events. f-branch statistic (f[b](C)) was calculated using Fbranch command, and its results were plotted using dtools.py. Although the default parameters for the f-branch statistic in Dsuite only consider fb with p < 0.01, we also assessed statistical significance of fb using a block Jack-knife approach by including the -Z parameter when running the f-branch statistic in Dsuite. A Z score |Z | > 3 was considered as significant. We generated diploid consensus sequences using the ‘mpileup’ command from the SAMtools ([237]http://samtools.sourceforge.net/) package, applying the recommended PSMC settings (-C 50, -O, -D 2*reads_depth, -d 1/3*reads_depth), and subsequently created input files with the ‘fq2psmcfa’ tool. For PSMC^[238]25 analysis, we set T[max] = 20, n = 64 (“4 + 50*1 + 4 + 6”) following Groenen et al.^[239]4. For MSMC2^[240]26, split-time estimation was conducted using phased genomes, following SNP calling, filtering and phasing as previously described^[241]7. Two masks were applied: one from ‘bamCaller.py’ of the MSMC-tools package and the other using Heng Li’s SNPable mask ([242]http://lh3lh3.users.sourceforge.net/snpable.shtml), which excludes non-mappable and non-unique sequences. MSMC2 was executed with four haplotypes (two from each population), calculating the relative cross coalescence rate (RCCR) with the “--skipAmbiguous” option to exclude unphased genomic regions. Time segment parameters were set to 64 (“4 + 50*1 + 4 + 6”) with 20 iterations. The RCCR values, ranging from 0 to 1, indicate population history, with RCCR approaching 1 suggesting the two populations were once a single population, and RCCR = 0.5 marking the estimated time of divergence between the populations. Finally, we used momi2^[243]29, which fits models to the site-frequency spectrum, to further assess the underlying population history based on six populations (EUW, RUE, VWB, SCW, NCW and RUA), assuming a mutation rate of 3.6 × 10^–9 per site per generation^[244]7. We used momi2 to contrast the fit of six evolutionary models proposed based on the above analyses results, e.g., model with no migration and other models where we fixed the gene flow from EUW and RUA to RUE and varied direction other migration we focused (Supplementary Fig. [245]5). We referred to the Akaike information criterion (AIC) to select model with smallest AIC as the optimal one. Then, we fitted 20 independent runs with different starting parameters and kept the model with the biggest log-likelihood value for the optimal model. Finally, we conducted 100 bootstrap calculations for the estimation of parameter range. Genome scan of selective sweeps To get more power to detect robust selective signal, we combined the SCW and VWB into a single gene pool representing South and East Asian wild boars (ASE) of comparative warm region, which is in contrast with RUA in the cold. The Fst statistics^[246]103 and π-ratio^[247]104 were used to identify signature of positive selection in RUA via VCFtools^[248]97. A third approach [the updated cross-population composite likelihood ratio test (XP-CLR), [249]https://github.com/hardingnj/xpclr]^[250]105 was performed for the comparison of RUA versus ASE to detect the selective signatures based on phased genotypes that were inferred by SHAPEIT4^[251]106. The command line was xpclr --input input.vcf --samplesA RUA_list --samplesB ASE_list –chr chrN --rrate 9.5e-9 --phased --maxsnps 500 –out outputFile. The recombination rate (--rrate 9.5e-9) of pig above was referred to Johnsson, et al. ^[252]107. In order to reliably identify signatures of selection and to discern selective sweeps from potential background divergence caused by bottleneck effects and from the comprehensive evidence of selection, we computed a single unbiased Composite Selection Score (CSS) by combining F[st,] π-ratio and XP-CLR estimates, following the procedure proposed by Randhawa, Khatkar, Thomson and Raadsma^[253]30. To compute the CSS, all loci are ranked for each statistic. For each locus, the fractional rank positions for each calculation (e.g., Fst), respectively, are converted into z-statistics, the mean of which corresponds to a single joint rank for all statistics. The corresponding p-value is retrieved from a standard normal distribution. The CSS score is then taken as–log[10]p. Protein-coding genes in the outlier windows were annotated using ANNOVAR software^[254]108. The specific mutation favored by selection in selective sweeps (identified by the CSS approach) was captured with phased genotypes using iSAFE (v1.1.1)^[255]34 with parameter --sample-cont ASE_list to set ASE as the control population and enabling the IgnoreGaps flag. In this analysis, we regarded the genotype of Sus verrucosus as the ancestral allele state because its high-quality sequencing. Transcriptomic analysis We collected 10 tissues (hypothalamus, liver, heart, muscle, skin, kidney, spleen, stomach, bone and subcutaneous fat) from 6 adult RUA individuals and 6 adult SCW individuals in winter, respectively. The sampling site for RUA was located in the Ternei District, Amgu village, in Primorsky Krai, Russia, while SCW was sampled in Jiangxi Province. Sampling was conducted in late December, mid-winter for both locations, with distinct climatic conditions: temperatures ranged from -14 °C to -7 °C in Primorsky Krai and from 5 to 12 °C in Jiangxi Province. There are some tissues of individuals excluded for transcriptomic analysis because of sample degradation (see Supplementary Data [256]3 for details). The age of all wild boars ranged approximately from 1 to 3 years, and detailed records of their sex are provided in Supplementary Data [257]3. RNA-seq reads (Illumina) were sequenced on the Illumina HiSeq 2500 platform. After quality control, the program HISAT2 (v2.1.0)^[258]109 was used to align sequencing reads to the reference genome. Annotated reference gene abundances were quantified using StringTie^[259]110. DESeq2^[260]111 was used to detect differential gene expression among different tissue comparisons between RUA and SCW with the criteria of padj < 0.05 and |log2FoldChange | > 1.2. Genes that were up-regulated in different tissues of RUA were used for GO and KEGG pathway enrichment analysis with R package clusterProfilter 4.0^[261]112. Estimating mitochondrial DNA (mtDNA) copy number from sequencing data First, we used SAMtools to obtain the coverage of each base in the genome from the aligned bam files. The average coverages for autosomal DNA and mtDNA were then calculated accordingly. In accordance with the analysis process of Yang et al. ^[262]113, we calculated the mtDNA copy number in various pig populations. Statistics and reproducibility Detailed explanations of the statistical methods used in this study can be found in the relevant sections of the Methods. For the selection scan, 42 RUA individuals and 25 ASE individuals were used. In the transcriptomic analysis, 6 individuals each of Southern wild boar and Russian wild boar were collected. Differentially expressed genes were identified using DESeq2, where the padj values were adjusted using the Benjamini-Hochberg method. Genes were considered significantly differentially expressed if padj < 0.05 and |log2FoldChange | > 1.2. When performing enrichment analysis using clusterProfiler 4.0, the Benjamini-Hochberg method is used to adjust the p-values. For figures, each data point corresponds to an individual animal. Ethics All procedures involving pigs were based on the care and use guidelines of experimental animals established by the Ministry of Science and Technology of China. The ethics committee of Jiangxi Agricultural University approved this study (JXAULL-2022-06-10). We have complied with all relevant ethical regulations for animal use. Reporting summary Further information on research design is available in the [263]Nature Portfolio Reporting Summary linked to this article. Supplementary information [264]Peer Review file^ (2MB, pdf) [265]Supplementary Information^ (11.9MB, pdf) [266]42003_2025_7536_MOESM3_ESM.pdf^ (38.7KB, pdf) Description of Additional Supplementary Files [267]Supplementary Data 1-24^ (5.8MB, xlsx) [268]Supplementary Data 25^ (3MB, xlsx) [269]Reporting summary^ (1.3MB, pdf) Acknowledgements