Investigation of single nucleotide polymorphisms in differentially expressed genes and proteins reveals the genetic basis of skeletal muscle growth differences between Tibetan and Large White pigs

Data utilization

For this study, we used whole-genome SNPs that displayed ΔAF between the TIBP and LWP breeds. These SNPs were derived from our previous study, which involved whole-genome sequencing of 44 TIBP and 29 LWP individuals. This dataset produced 16.47 billion raw reads, with a mean depth of 9.13× per individual and average genome coverage of 98.96%. A total of 21,594,848 and 15,210,134 SNPs were identified for TIBP and LWP, respectively. In total, 21,767,938 SNPs were obtained from the 73 individuals for subsequent analysis. SNPs were categorized based on their specificity and predominance in each breed. SNPs with allele frequencies of 0.5 or higher in TIBP and nonexistent in LWP were classified as TIBP-specific. Those with allele frequencies greater than 0.5 in TIBP and less than 0.5 in LWP, accompanied by an absolute ΔAF equal to or exceeding 0.5, were termed TIBP-predominant SNPs. Conversely, SNPs with allele frequencies of 0.5 or higher in LWP and nonexistent in TIBP were designated LWP-specific. SNPs with allele frequencies greater than 0.5 in LWP, less than 0.5 in TIBP, and ΔAF between LWP and TIBP of 0.5 or more were labeled as LWP-predominant SNPs.

Shang et al [16] detected the transcriptome and preteome profiles in longissimus dorsi muscle tissues of nine embryos from two pregnant TIBP and LWBP sows at 60 days after insemination. They identified 3,858 DEGs and 830 DEPs, respectively, between the TIBP and LWP breeds. In this study, we integrated the identified 3,858 DEGs and 830 DEPs into our analysis. This integration aimed to enrich our understanding of the genetic variations of genes involved in skeletal muscle development divergenece between these two pig breeds.

Identification of TIBP- and LWP-specific and predominant SNPs in DEGs and DEPs

In our study, we meticulously annotated specific and predominant SNPs for TIBP and LWP breeds and assigned them to corresponding genes, referred to as “candidate genes.” We then conducted a comprehensive comparison of the candidate genes, DEGs, and DEPs. This process enabled us to identify genes that overlapped between the candidate genes and DEGs and DEPs. To visualize the overlap and relationships among these genes, we used OmicStudio tools (available at https://www.omicstudio.cn/tool), creating a Venn diagram for an intuitive representation of our findings. Additionally, we used custom-developed Python scripts to extract SNPs within overlapping DEGs and DEPs.

Functional enrichment analysis of overlapping DEGs

To delve deeper into the functional implications of the overlapping DEGs, we conducted a gene ontology (GO) enrichment analysis using g:Profiler tools (available at https://biit.cs.ut.ee/gprofiler/orth) [17]. In this analysis, we set organism to "Sus scrofa", statistical domain scope set to "only annotated genes", significance threshold set to "Benjamini-Hochberg FDR", and the user threshold set to "0.05".

Subsequently, we performed Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis using the DAVID database (https://david.ncifcrf.gov/tools.jsp) [18], a widely recognized platform for high-throughput gene function analysis. Again, the organism was specified as "Sus scrofa," and we adhered to the default settings for the other parameters. This approach enabled us to identify potential genes involved in skeletal muscle development, as indicated by their enrichment in the relevant GO terms and KEGG pathways.

We investigated the allele frequencies of SNPs within these identified genes to further refine our insights. This was achieved using in-house custom-developed scripts, which allowed us to identify specific genetic variations that potentially drive skeletal muscle development.

Functional enrichment analysis of overlapping DEPs

We conducted a GO enrichment analysis to elucidate the functional roles of the overlapping DEPs in skeletal muscle growth. This was achieved using the g:Profiler tool (available at https://biit.cs.ut.ee/gprofiler/orth) [17]; the organism was set to "Sus scrofa", the statistical domain scope was set to "only annotated genes", the significance threshold was set to "Benjamini-Hochberg FDR", and the user threshold was set to "0.05".

Furthermore, we performed a KEGG pathway analysis using the DAVID database (https://david.ncifcrf.gov/tools.jsp) [18]. For this analysis, the organism was again set to "Sus scrofa", and the default settings were maintained for other parameters. This approach allowed us to identify DEPs that were significantly enriched in GO terms and KEGG pathways pertinent to skeletal muscle development, suggesting their potential role in this process.

We investigated the allele frequencies of SNPs contained within them to deepen our understanding of these DEPs. This task was accomplished using custom-developed scripts, enabling a detailed examination of genetic variations that potentially influence skeletal muscle development.

TIBP- and LWP-specific and predominant SNPs in DEGs and DEPs

In our previouly study, 2,893,106 and 813,310 specific and predominant SNPs in the TIBP and LWP populations were identified and annotated to 24,560 genes. In this study, we referred to these annotated genes as “candidate genes”. Shang et al [16] reported a significant number of DEGs and DEPs, specifically 3,858 DEGs and 830 DEPs, in the longissimus dorsi muscle tissues from the TIBP and LWP breeds at 60 days of embryonic development. Upon comparing these candidate genes with DEGs and DEPs, we discovered that 3,499 DEGs and 628 DEPs overlapped with candidate genes (Figure 1A). This overlap represented 91.27% of the total DEGs and 75.67% of the total DEPs, indicating that a major part of these genes and proteins contained TIBP- and LWP-specific and predominant SNPs. Specifically, 3,499 DEGs contained 924,051 specific and predominant SNPs, including 270,189 TIBP-specific and 449,608 predominant SNPs, along with 190,210 LWP-specific and 14,045 predominant SNPs. Furthermore, 628 DEPs contained 109,625 specific and predominant SNPs, comprising 32,938 TIBP-specific and 53,692 predominant SNPs, in addition to 1,983 LWP-specific and 21,014 predominant SNPs. We observed that these specific and predominant SNPs were distributed across all chromosomes (Figure 1B and 1C). The ΔAF values of these SNPs ranged from 0.5 to 1.0, with the majority falling within the 0.5 to 0.6 range (illustrated in Figure 1D–1G). Furthermore, we identified 263 genes shared among the DEGs, DEPs, and candidate genes (Figure 1A).

Figure 1

Summary of TIBP- and LWP-specific and predominant SNPs. (A) A Venn diagram illustrating the overlap of genes among candidate genes, DEGs, and DEPs; (B) and (C) the chromosomal distribution of TIBP- and LWP-specific and predominant SNPs found within DEGs and DEPs; and (D)–(G) the counts of SNPs with ΔAF ranging from 0.5 to 1.0 for TIBP- and LWP-specific and predominant SNPs harbored by DEGs and DEPs. TIBP, Tibetan pig; LWP, Large White pig; SNPs, single nucleotide polymorphism; DEGs, differentially expressed genes; DEPs, differentially expressed proteins; AF, allele frequency differences.

Investigation of SNPs in DEGs related to skeletal muscle development

We conducted a GO enrichment analysis of the overlapping DEGs to investigate SNPs present in the DEGs associated with skeletal muscle growth. This analysis revealed that 3,499 overlapping genes were significantly enriched in 1,481 biological processes, 246 cell components, and 151 molecular functions (Supplementary Table S1). Focusing specifically on skeletal muscle development, a type of striated muscle, we focused on the biological processes relevant to this area. We identified seven striated muscle biological process-related GO terms, encompassing striated muscle cell differentiation, development, tissue development, contraction, adaptation, hypertrophy, and negative regulation of muscle contraction. Among these categories, 145 genes were identified, including a significant number of SNPs: 5,983 TIBP-specific and 16,894 predominant SNPs, along with 433 LWP-specific and 7,072 predominant SNPs (Supplementary Table S2). These SNPs were classified into 13 distinct types, including missense, synonymous, 5′ prime untranslated regions (UTR), and 3′ prime UTR variations. The majority of these SNP were variations in the intronic region, accounting for 66.2% of the total. Supplementary Table S3 presents a detailed breakdown of each SNP type, providing a thorough understanding of the SNP landscape in DEGs associated with skeletal muscle development.

We identified 145 genes, including two members of the MRF family (MYF5 and MYF6) and three main members of the MEF family (MEF2A, MEF2C, and MEF2D), which are key to skeletal muscle development. Additionally, genes such as MSTN, which act as negative regulators of muscle development, were identified. SNPs located in the exonic, splice, and regulatory regions, including synonymous, missense, splice region, 3′ prime UTR, and 5′ prime UTR, as well as downstream and upstream variations, can significantly influence gene structure and function. Breed-specific SNPs are crucial in defining breed-specific traits. Therefore, we focused on genes containing TIBP-specific SNPs in regulatory and/or coding regions, especially those with allele frequencies exceeding 0.80. This led to the discovery of several genes harboring multiple regulatory and/or coding region TIBP-specific SNPs with allele frequencies greater than 0.8. Table 1 presents detailed information on these SNPs. The genes identified include MYF5, MYOF, and ASB2 (implicated in striated muscle cell differentiation); PDE9A (associated with striated muscle hypertrophy); SDC1, PDGFRA, MYOM2, ACVR1, and ZIC3 (involved in striated muscle cell development); COL11A1, TGFBR1, EDNRA, and TGFB2 (contributing to striated muscle tissue development); PDE4D, PGAM2, GRK2, and SCN4B (related to striated muscle contraction); and CACNA1S and MYL4 (associated with striated muscle adaptation).

Table 1

Tibetan pigs-specific single nucleotide polymorphisms in regulatory and/or coding regions harbored by differentially expressed genes related to skeletal muscle development

KEGG pathway analysis revealed significant enrichment of over-represented genes across 118 distinct pathways, as detailed in Supplementary Table S4. These pathways encompass critical biological processes, including energy production (such as metabolic pathway, carbon metabolism, citrate cycle, and pyruvate metabolism), fat metabolism (fatty acid metabolism and fatty acid degradation), adipocyte differentiation (PPAR signaling pathway), actin dynamics (regulation of actin cytoskeleton), and amino acid metabolism (biosynthesis of amino acids, valine, leucine, and isoleucine degradation, arginine biosynthesis, arginine, and proline metabolism). Several signaling pathways, including insulin [11], AMPK [19], FoxO [10], PI3K-Akt [20], ERBB [21], HIF-1 [22], MAPK [23], and P53 signaling pathways [24], play roles in various biological processes, notably in protein synthesis and degradation via regulation of the mTOR signaling pathway. Consequently, we identified 22 genes involved in the mTOR signaling pathway. This set included genes such as IGF1, INSR, IRS3, PIK3R2, SGK1, AKT2, SKP2, MDM2, RPS6KB2, RHEB, EIF4EBP1, SHC1, RAF1, FOXO1, MAP2K2, MKNK2, STK11, AKT1S1, RB1CC1, ULK1, ERBB2, and GAB1, which collectively harbor 3,915 specific and predominant SNPs (Supplementary Table S5). Supplementary Table S6 presents detailed classifications of these SNP types. Further analysis of allele frequencies revealed that genes such as IGF1, FOXO1, and EIF4EBP1 possess functionally significant TIBP-specific SNPs with allele frequencies exceeding 0.80, as shown in Table 2. Notably, IGF1 includes several SNP types in regulatory regions, including the 5′ prime UTR, upstream, 3′ prime UTR, and downstream regions.

Table 2

Tibetan pig-specific single nucleotide polymorphisms in the regulatory region harbored by differentially expressed genes involved in the mTOR signaling pathway

Investigation of SNPs in DEPs related to skeletal muscle development

Our GO enrichment analysis of the 628 overlapping DEPs identified significant enrichment across 270 biological processes, 177 cell components, and 95 molecular function terms, as detailed in Supplementary Table S7. This analysis revealed four biological processes related to the GO terms associated with striated muscle, including striated muscle cell development, hypertrophy, contraction, and adaptation. Within these categories, 23 genes were identified. These genes collectively harbored 3,114 specific and predominant SNPs, which are cataloged in Supplementary Table S8. These SNPs spanned 12 distinct SNP types, with a notable majority being intergenic region variations, accounting for 57.1%, as outlined in Supplementary Table S9. An examination of allele frequencies highlighted several significant findings: one upstream TIBP-specific SNP in the MYL4 genes with 0.93 allele frequency, a splice region variant and intron TIBP-specific SNP in the MYOM2 gene with an allele frequency of 0.86, and two 5′ prime UTR plus one 5′ prime UTR premature start codon gain TIBP-specific SNPs in PGAM2, averaging an allele frequency of 0.84.

KEGG pathway analysis revealed significant enrichment in 51 over-represented pathways, as detailed in Supplementary Table S10. These pathways encompass a broad range of biological functions, including energy production (e.g., metabolic pathways, carbon metabolism, and citrate cycle), fat metabolism (such as fatty acid metabolism and fatty acid degradation), and amino acid metabolism (including biosynthesis of amino acids, valine, leucine, and isoleucine degradation, arginine biosynthesis, arginine, and proline metabolism). Furthermore, pathways integral to protein synthesis and degradation were highlighted, notably the insulin and HIF-1 signaling pathways. Within the realm of these pathways, genes such as EIF4EBP1, EIF4E2, and EIF4E have been identified as being involved in the mTOR signaling pathway.

SNP analysis of shared genes in DEGs and DEPs linked to skeletal muscle development

In our study, we delved deeper into the shared genes found between DEGs and DEPs that play roles in biological processes associated with striated muscle development. The results revealed that 18 genes were overlapped. These are CASQ1, MYOM2, TPM1, ACTA1, CAPN3, NEB, and MYOZ1 (involved in striated muscle cell development); TCAP, RAP1GDS1, ATP2B4, LMCD1, and TWF1 (associated with striated muscle hypertrophy); and TNNC2, MYL4, ATP2A1, PGAM2, HOMER1, and STAC3 (related to striated muscle contraction). A total of 2,930 specific and pedominant SNPs were identified within these genes (Supplementary Table S11). Notably, several of these genes harbored functionally significant specific and predominant SNPs, with ΔAF exceeding 0.80. These included genes such as TPM1, CAPN3, RAP1GDS1, MYOZ1, MYOM2, and PGAM2, as detailed in Table 3. Specifically, the CAPN3 gene contained four synonymous and three 5′ prime UTR TIBP-predominant SNPs; the MYOM2 gene had five synonymous and one splice region variant and intron TIBP-predominant SNPs; and the PGAM2 gene included one 5′ prime UTR premature start codon gain, two 5′ prime UTR TIBP-specific SNPs, four 5′ prime UTR SNPs, one 5′ prime UTR premature start codon gain and one synonymous LWP-predominant SNP.

Table 3

Specific and predominant single nucleotide polymorphisms in regulatory and/or coding regions harbored by shared genes among DEGs and DEPs related to skeletal muscle development

Additionally, EIF4EBP1 emerged as a shared gene among both DEGs and DEPs, as identified in our KEGG pathway analysis. The EIF4EBP1 gene was found to contain 22 TIBP-specific SNPs, comprising 20 intronic variations and two upstream variations. Within the TIBP population, the SNPs exhibited an average allele frequency of 0.76, highlighting their genetic distinctiveness and potential functional impact in this breed.

Potentially crucial DEGs related to skeletal muscle growth

The MRF and MEF families are pivotal in muscle cell differentiation and development and play central roles in the regulation of myogenesis [9]. Notably, in our study, two members of the MRF family, MYF5 and MYF6, and three principal members of the MEF family, MEF2A, MEF2C, and MEF2D, emerged as DEGs. Intriguingly, these genes harbored multiple TIBP- and LWP-specific and predominant SNPs. Moreover, MSTN, a member of the transforming growth factor beta (TGF-β) superfamily produced by skeletal muscle, acts as a key inhibitor of muscle growth. The occurrence of myostatin mutations across various mammalian species, which lead to muscle hypertrophy, underscores its vital role in modulating muscle development and size. Our study further identified MSTN as a DEG containing several TIBP- and LWP-specific and predominant SNPs. This finding aligns with previous research, underscoring the significance of these genes in muscle fiber formation and potentially elucidating the observed disparities in muscle growth rates between the TIBP and LWP breeds.

The discovery of SNPs within the regulatory and/or coding regions of DEGs and DEPs associated with muscle growth is a key achievement of our study. Particularly noteworthy were the TIBP- and LWP-specific and predominant SNPs that exhibited marked ΔAF between the TIBP and LWP breeds. These SNPs potentially contribute to the slower muscle growth observed in TIBP than in LWP. SNPs located upstream, 3′ prime UTR, 5′ prime UTR, downstream, and splice regions have the potential to modulate gene expression levels and patterns [25]. Their influences on post-transcriptional regulation, particularly affecting mRNA stability and translation efficiency, may play a crucial role in driving phenotypic disparities in muscle development observed between the TIBP and LWP breeds. Furthermore, SNPs in coding regions, such as missense variations, can alter amino acid sequences [26], thereby influencing the structure and functionality of proteins vital for muscle development. Even synonymous SNPs, traditionally thought to be silent, can affect gene expression and protein folding, as indicated in a recent study [27].

In this study, we observed that GRK2, ASB2, PDGFRA, ACVR1, PGAM2, IGF1, and FOXO1 were DEGs containing several TIBP-specific SNPs in regulatory and/or coding regions with allele frequencies greater than 0.8. GRK2 plays a crucial role in skeletal muscle development and differentiation, primarily by regulating key signaling pathways, such as p38MAPK and Akt. Studies in various models, including Drosophila and mice, have shown that alterations in GRK2 levels can lead to defects in muscle differentiation, changes in muscle fiber size, and impaired muscle function. Thus, balanced GRK2 activity is essential for proper myogenic processes and skeletal muscle growth [28]. ASB2 functions as a negative regulator of skeletal muscle mass, influenced by the TGF-β signaling network, and follistatin-based interventions, which modulate TGF-β network activity, lead to decreased ASB2 expression, promoting muscle hypertrophy [29].

PDGFRA plays a pivotal role in directing stem cell fate toward the paraxial mesoderm, which is essential for the formation of myogenic precursor cells [30,31]. ACVR1 plays a significant role in myogenic differentiation, as demonstrated by its mutation (Acvr1R206H/+) negatively affecting this process. The mutation in fibro/adipogenic progenitors (FAPs) impairs their ability to support the myogenic potential of muscle satellite cells (MuSCs), leading to reduced fusion and formation of myofibers. Conversely, healthy FAPs can ameliorate the impaired myogenic morphology of mutant MuSCs, and the ACVR1 signaling pathway, especially in the context of the transforming growth factor beta/bone morphogenetic protein (TGFβ/BMP) pathway, is crucial, as the dysregulation of components such as BMP2 and Noggin in FAP-conditioned media is associated with impaired muscle regeneration and differentiation [32]. PGAM2, a key glycolytic enzyme highly expressed in skeletal muscle, is critical for myogenic differentiation. Mutations in PGAM2, such as the K176R alteration, lead to impaired myogenic differentiation and reduced cellular glycolysis and mitochondrial function, as evidenced in CRISPR-engineered C2C12 myogenic cells [33].

IGF1 plays a pivotal role in regulating the growth and maintenance of skeletal muscles. It promotes muscle protein synthesis through the PI3K/Akt/mTOR and PI3K/Akt/GSK3β pathways and inhibits muscle atrophy by suppressing E3 ubiquitin ligases via the PI3K/Akt pathway, thereby reducing protein degradation in the ubiquitin-proteasome system [20]. FoxO1, a member of the FoxO family, is a crucial transcription factor in skeletal muscle that plays a dual role in muscle growth and differentiation. While some studies suggest that FoxO1 acts as a negative regulator of skeletal muscle differentiation, others highlight its importance in myoblast fusion [34]. Furthermore, the co-occurrence of SNPs in both the regulatory and coding regions underscores the genetic complexity of muscle development in pigs. This complexity suggests a multilayered regulation of key genes that influence muscle growth through both transcriptional and post-transcriptional mechanisms.

Shared genes between DEGs and DEPs potentially crucial for skeletal muscle growth

In addition, our focus was extended to genes that were both DEGs and DEPs related to skeletal muscle growth, and possessed TIBP- and LWP-specific and predominant SNPs in regulatory and/or coding regions with ΔAF greater than 0.8. Concurrent changes in the levels of the gene and its protein indicate that the function of this gene is critical under certain conditions [35]. In our study, CAPN3, MYOM2, and PGAM2 were identified in this category.

CAPN3 is the only muscle-specific calpain that plays an important role in the promotion of calcium release from skeletal muscle fibers, calcium uptake by the sarcoplasmic reticulum, muscle formation, and muscle remodeling [36]. We discovered four synonymous, three 5′ prime UTR, and three upstream TIBP-predominant SNPs in the CAPN3 gene, suggesting a significant role in the growth disparities between the TIBP and LWP breeds.

Five synonymous TIBP-predominant SNPs were detected in MYOM2 in this study. MYOM2 plays a significant role in skeletal muscle function. It is involved in the structure and integrity of the muscle sarcomere, particularly at the M-line, where it contributes to the assembly and stabilization of thick filaments. MYOM2 interacts with other sarcomeric proteins to maintain the proper alignment and function of the sarcomere during muscle contraction. Mutations or dysregulation of MYOM2 can affect muscle strength and may be implicated in various muscle disorders [37].

PGAM2, a crucial glycolytic enzyme highly expressed in skeletal muscle and essential for myogenic differentiation [33], exhibited two 5′ prime UTR premature start codon gains, six 5′ prime UTR, and one synonymous TIBP-specific and LWP-predominant SNPs. These findings underscore the potential critical role of these genes in the differences in primary fiber development between the TIBP and LWP breeds, highlighting their importance in skeletal muscle growth and development.

Toward a holistic understanding of genetic determinants in skeletal muscle development across life stages in TIBP and LWP breeds

This study synergistically combined DEGs and DEPs from the critical stage of primary myofiber formation with genomic data to uncover potential causative genes for the differences in skeletal muscle growth between the TIBP and LWP breeds. However, it is important to recognize that skeletal muscle development is a polygenic trait exhibiting varied characteristics across different stages of life. For example, during the prenatal phase, the development of skeletal muscle primarily involves an increase in myofiber numbers, beginning with the formation of primary myofibers, followed by secondary myofibers [5]. Conversely, postnatal muscle growth is predominantly characterized by hypertrophy of existing fibers [38,39]. Therefore, a joint analysis that encompasses DEGs and DEPs from multiple developmental stages, combined with genomic data, would enable more thorough and insightful exploration. This integrative approach is crucial for identifying the genes responsible for the observed differences in skeletal muscle development between the TIBP and LWP breeds, thereby providing a complete understanding of the genetic factors influencing muscle growth across different life stages.