Genome-Wide Selection Signals of the Min Pig Indicate That Genes Related to Fat and Energy Metabolism Play a Crucial Role in Responding to Cold Stress

1. Introduction

Pork, as one of the most consumed meats globally, is a significant source of protein and nutrients in the human diet [1]. Indigenous pig breeds hold unique value in the pig industry and breeding, as they have undergone long-term natural selection and adaptive evolution, enabling efficient production under specific environmental conditions and providing valuable genetic resources for the sustainable development of the global pig industry. The Min Pig (MZ), a typical representative of Chinese indigenous pig breeds, has developed distinct genetic traits through prolonged natural selection and domestication, particularly excelling in cold adaptability [2].

With the completion of the Swine Genome Project, breeders worldwide have successfully identified multiple key genes influencing economic traits in pigs using high-density genotyping and resequencing technologies. For instance, genes related to growth rate (e.g., FASN, ELOVL6, SSC7, SSC2) [3,4,5], lean meat percentage (e.g., NCK2, ACSL1, PPARG) [6,7,8] , meat quality (e.g., GHRHR, SSC18, CDYL2) [7,9,10] , and litter size (e.g., IGFBP2, UBE3A, FAM135B) have been pinpointed [11,12,13].

In recent years, a series of candidate genes related to environmental adaptability in pigs have been identified. For example, FBN1 was confirmed that closely related to skin thickness in pigs and may involve in their environmental adaptability [14,15]. Previous studies identified several candidate genes of physiological responses of pigs under heat stress, including GHR, TEAD4, NNT, ERBB4, FKBP1B, and NFATC2 [16]. Additionally, current study indicated that the haplotypes of the BTF3 gene significantly differ between Chinese and Western pig breeds, influencing fat tissue formation by promoting adipocyte proliferation and regulating intramuscular fat deposition [17]. Furthermore, the upregulation of glucose metabolism-related genes (e.g., ACSS1, HK3, PGM2, PCK1) under cold stress environments was confirmed to enhance glycolysis, gluconeogenesis, and glycogenolysis, increasing energy production to counteract cold exposure [18]. In this study, a genome-wide selective sweep analysis (GWSA) of MZ was performed to systematically identify selected genes relative to their cold adaptability, which will be helpful to deepen the genetic basis understanding of cold adaptability in pigs.

2. Materials and Methods

In this study, total of 30 MZ blood samples, with no familial relationships within three generavation farm in Lanxi City, Heilongjiang Province of China (E 126°16′, N 46°15′). Genomic DNA was extracted using TaKaRa MiniBEST Universal Genomic DNA Extraction Kit Ver.5.0 (Takara, Japan). Sequencing libraries were constructed by Annoroad® Universal DNA Library Preparation Kit v2.0 (Illumina®, San Diego, CA, USA) and sequenced on the BGISEQ-500 platform (Beijing Genomics Institute, China) with a sequencing depth of ~10× per individual. Additionally, the publicly available genomic data of 42 samples from 7 breeds without cold adaptability were obtained from the Sequence Read Archive (SRA) database (Table S1) as control group for subsequent GWSA. These populations included 3 indigenous pigs in southern China, such as Chenghua (CH, n=5), Rongchang (RC, n=7) and Meishan (MS, n=5) and 4 commercial breeds, such as Duroc (DL, n=6), Yorkshire (YK, n=7), Landrace (CB, n=5), and French Large White (DB, n=7).

Integration with publicly available data and 30 MZ genome re-sequence data produced from this study, a total of 1566.8 GB of raw sequencing data were obtained. The raw data were quality-controlled and filtered using Fastp (v0.23.4) with default parameters to obtain high-quality clean reads (HQRs). Subsequently, all HQRs were aligned to the pig reference genome (GCF_000003025.6, Sscrofa 11.1) using bwa-mem2 (v2.2.1) for further analysis. Single nucleotide polymorphisms (SNPs) were identified using the HaplotypeCaller module in the Genome Analysis Toolkit (GATK, version 4.5.0.0), and variant quality filtering was performed using the VariantFiltration module (McKenna et al. 2010). The filtering criteria were set as follows: quality depth (QD) <2.0, variant quality (QUAL) <30.0, strand bias ratio (SOR) >3.0, Fisher's exact test (FS) >60.0, mapping quality (MQ) <40.0, mapping quality rank sum test (MQRankSum) <-12.5, and read position rank sum test (ReadPosRankSum) <-8.0. Further filtering of autosomal SNPs was performed using vcftools (v0.1.16) [19] with the following parameters: only biallelic variants retained (--min-alleles 2 --max-alleles 2), genotype missing rate <0.05 (--max-missing 0.95), and minimum allele frequency > 0.05, (--maf 0.05). Finally, candidate SNPs were functionally annotated using Annovar software [20].

The GWSA was conducted on MZ (case) and other non-cold-adapted breeds (control) using 3 statistical algorithms: fixation index (FST) [21], nucleotide diversity difference (θπ ratio) [22] , and cross-population composite likelihood ratio test (XP-CLR) [23]. The analysis of these parameters was estimated with a sliding window of 40 kb and a step size of 20 kb. The top 5% of windows (thresholds set as follows: FST≥ 0.412596, θπ ratio≥ 2.125852484, and XP-CLR≥ 22.08201682) were selected from each parameter, and regions overlapping across all three methods were considered candidate selective regions (CSRs). The candidate genes (CDGs) were defined as the genes overlapped in above CSRs. Subsequently, the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed on CDGs using online tool KOBAS (http://kobas.cbi.pku.edu.cn/), with a significance threshold set at a corrected P-value <0.05.

3. Results and Discussion

A total of 1,566.8 GB of raw sequencing data were obtained. Across all 72 individuals, 60,652,730 SNPs were identified, and 23,991,936 high-quality SNPs were retained for subsequent analysis after quality control. The selection signal results indicated that the windows with the highest FST values were located on chromosomes CHR4 and CHR5 (Figure 1A), while the windows with the greatest θπ ratio differences were found on CHR4 and CHR8 (Figure 1B). The windows with the highest XP-CLR values were located on CHR4 and CHR15 (Figure 1C). Based on the intersection of top 5% windows from these three parameters, total 58,771 SNPs were obtained from 357 intersected windows, which located in exonic, intronic, 3’UTR, 5’UTR, and upstream/downstream regions. Finally, total 226 candidate genes were identified originating from 357 intersecting windows.

Additionally, gene functional annotation results revealed that 68 of 226 CDGs were enriched in 183 KEGG pathways, with 24 KEGG pathways significantly enriched (corrected-P < 0.05, Table 1), including the cGMP-PKG signaling pathway (NFATC2, MAP2K2, MYLK4, CREB3L3, AKT3, PRKG1), PI3K-Akt signaling pathway (ERBB4, CREB3L3, MAP2K2, GNG7, AKT3, MDM2, ANGPT1, BDNF), and Kaposi sarcoma-associated herpesvirus infection (NFATC2, MAP2K2, PREX1, GNG7, AKT3, LYN). Furthermore, numerous CDGs (e.g., AKT3, PRKG1, CREB3L3, ACSF3, and NDUFS7) were enriched in pathways related to fat synthesis and energy metabolism, such as Fatty acid biosynthesis, Carbohydrate digestion and absorption, Adipocytokine signaling pathway, Regulation of lipolysis in adipocytes, Glucagon signaling pathway, Insulin resistance, and Thermogenesis, etc.

As well-known, fat plays a crucial role in resisting low-temperature environments, serving as a key substance for providing energy, maintaining body temperature, and enhancing cold stress resistance. This study identified a large number of genes related to fat metabolism and synthesis that are under selection in MZ. For example, ACSF3 is a malonyl-CoA synthetase involved in the synthesis of fatty acids in mammalian mitochondria [24], and its activity also affects the malonylation level of ACOT7, thereby regulating fat content in animal muscles [25,26]. ATP6V1H is the H subunit of the V1 subunit of V-ATPase, participating in the assembly and functional regulation of V-ATPase, which is crucial for maintaining cellular acid-base balance and various physiological processes [27]. Research has found that the absence of ATP6V1H increases endoplasmic reticulum stress in pancreatic β-cells and exacerbates glucose tolerance impairment caused by fatty acids [28]. This indicates that ATP6V1H plays a key role in cellular metabolism and stress response, and its dysfunction may affect the body's adaptation to cold stress. NMRK2 is a muscle-specific β1 integrin-binding protein that activates the SIRT1 pathway by increasing NAD+ levels, thereby promoting fatty acid oxidative metabolism and enhancing cold resistance in animals [29,30], and it is potentially related to the development and regulation of goat adipose tissue [31]. The PIP5K1C gene encodes phosphatidylinositol-4-phosphate 5-kinase, which is involved in phospholipid metabolism and promotes the synthesis of diacylglycerol (DAG) and phosphatidylinositol (PI) [32]. The activity of PIP5K1C affects cell signaling pathways, especially those related to cell migration and signaling [33,34]. Previous studies have confirmed that the absence of PIP5K1C in adipocytes reduces the phosphorylation levels of PI3K and AKT, thereby inhibiting the activity of the PI3K/AKT signaling pathway and promoting cell apoptosis [35] .

Moreover, GAMT, as a crucial enzyme in the synthesis of creatine, regulates the creatine-phosphocreatine system to maintain high-energy phosphate reserves, providing rapid ATP supply to muscle and nerve tissues under low-temperature conditions [36] . Research indicates that GAMT deficiency decreases the efficiency of muscle energy metabolism, exacerbating metabolic imbalances caused by cold stress [37]. Its involvement in creatine metabolism indirectly enhances thermogenic capacity by modulating mitochondrial function, thereby improving cold adaptation [38] . The AMPD3 activates the AMPK pathway by regulating the AMP/ATP ratio, promoting fatty acid oxidation and thermogenic metabolism[39,40]. Overexpression of AMPD3 was also reported to increase mitochondrial oxygen consumption rates in muscles[41].Studies suggest that overexpression of CREB3L3 enhances lipolysis, ketogenesis, and insulin sensitivity, increasing energy expenditure and improving metabolic indicators, revealing its critical regulatory role in energy metabolism and thermogenic processes[42,43,44].

PRKG1, as a cGMP-dependent protein kinase, was suggested to play a critical role in the vasoconstrictive response induced by cold exposure. For example, study indicates that PRKG1 regulates the localization and function of the Rho A protein through its phosphorylation, to decreasing the phosphorylation of myosin light chains, thereby promoting vasodilation [45]. Additionally, PRKG1 facilitates the hydrolysis of triglycerides to release free fatty acids and glycerol[46]. Calcium-dependent neutral proteases, such as Calpain-1, also play a crucial role in the regulation of mitochondrial function and cellular stress responses, with their activity being associated with mitochondrial calcium overload. Inhibiting Calpain-1 has been shown to enhance mitochondrial function and reduce cellular damage[47], thereby maintaining thermogenic efficiency. NDUFS7 is the core catalytic subunit of mitochondrial complex I (NADH dehydrogenase), directly involved in the electron transfer from NADH to ubiquinone [48]. It was reported that an inhibitor (DX2-201) targets NDUFS7 and inhibits mitochondrial function and oxidative phosphorylation by suppressing the activity of complex I [49]. Furthermore, the deficiency of NDUFS7 leads to a decreased cell proliferation, increased cell death, and increased susceptibility to oxidative stress [50]. It is worth noting that studies have shown that exposure of animals with constant temperature to low environmental temperatures is associated with oxidative stress in various body tissues [51].

ACOT8 (Acyl-CoA Thioesterase 8) is a member of the acyl-CoA thioesterase superfamily and acts as a peroxisomal acyl-CoA thioesterase capable of hydrolyzing various acyl-CoA substrates[52]. Study indicates that ACOT8 facilitates the conversion of acetyl-CoA into acetate and CoA, thereby promoting the reutilization of CoA and supporting the ongoing process of fatty acid oxidation (FAO) during energy stress. Furthermore, ACOT8 is upregulated enhancing the function of HMGCS2 (3-Hydroxy-3-Methylglutaryl-CoA Synthase 2) to promote ketogenesis, thus supplying fuel to extrahepatic tissues and enhancing the body's cold tolerance[53]. Lmpad1 is a Golgi-resident protein that regulates Golgi morphology and vesicular transport[54]. Under low-temperature conditions, IMPAD1 may maintain Golgi function, facilitating the glycosylation of mitochondrial-related enzymes, which optimizes ATP production efficiency, decrease ROS levels and alleviates energy shortages [55].

B4GALT5 (β-1,4-galactosyltransferase 5) belongs to the β-4-galactosyltransferase gene family associated with the Golgi apparatus and plays a crucial role in the biosynthesis of glycosphingolipids and glycans. It primarily catalyzes the attachment of galactose to glycosphingolipid or glycoprotein receptors, such as GlcNAc, Glc, and Xyl, through β-1,4-glycosidic bonds[56]. Studies indicate that B4GALT5 is present not only on the cell surface but also within the Golgi complex, where it functions as an adhesion molecule involved in matrix interactions, cellular spreading and migration, and signal transduction cascades[57,58]. This suggests that B4GALT5 may enhance the cold stability of cell membranes and adhesion structures through glycosylation modifications, thereby promoting the cellular migratory capacity associated with cold injury repair and optimizing the efficiency of stress signaling pathways.

4. Conclusion

In this study, a series of genes involved in fat metabolism, energy homeostasis, and mitochondrial function that are under selection in MZ were confirmed by GWSA analysis. These findings not only deepen our understanding to genetic basis of cold adaptation, but also provide support for further development of environmental adaptation breeding markers for pigs.