Skip to main content

ORIGINAL RESEARCH article

Front. Genet., 11 April 2022
Sec. Livestock Genomics
This article is part of the Research Topic Marbling Fat in Livestock View all 13 articles

Integrative Analysis of Nanopore and Illumina Sequencing Reveals Alternative Splicing Complexity in Pig Longissimus Dorsi Muscle

Ze ShuZe ShuLigang Wang
Ligang Wang*Jinbu WangJinbu WangLongchao ZhangLongchao ZhangXinhua HouXinhua HouHua YanHua YanLixian WangLixian Wang
  • Key Laboratory of Farm Animal Genetic Resources and Germplasm Innovation of Ministry of Agriculture of China, Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, Beijing, China

Alternative splicing (AS) is a key step in the post-transcriptional regulation of gene expression that can affect intramuscular fat (IMF). In this study, longissimus dorsi muscles from 30 pigs in high- and low- IMF groups were used to perform Oxford Nanopore Technologies (ONT) full-length sequencing and Illumina strand-specific RNA-seq. A total of 43,688 full-length transcripts were identified, with 4,322 novel genes and 30,795 novel transcripts. Using AStalavista, a total of 14,728 AS events were detected in the longissimus dorsi muscle. About 17.79% of the genes produced splicing isoforms, in which exon skipping was the most frequent AS event. By analyzing the expression differences of mRNAs and splicing isoforms, we found that differentially expressed mRNAs with splicing isoforms could participate in skeletal muscle development and fatty acid metabolism, which might determine muscle-related traits. SERBP1, MYL1, TNNT3, and TNNT1 were identified with multiple splicing isoforms, with significant differences in expression. AS events occurring in IFI6 and GADD45G may cause significant differences in gene expression. Other AS events, such as ONT.15153.3, may regulate the function of ART1 by regulating the expression of different transcripts. Moreover, co-expression and protein-protein interaction (PPI) analysis indicated that several genes (MRPL27, AAR2, PYGM, PSMD4, SCNM1, and HNRNPDL) may be related to intramuscular fat. The splicing isoforms investigated in our research provide a reference for the study of alternative splicing regulation of intramuscular fat deposition.

Introduction

As pork is the main meat resources for humans, its palatability has always been an important concern. The content of intramuscular fat (IMF) in meat is strongly correlated with the overall meat quality. Increased IMF content not only promotes the sensory attributes (Alfaia et al., 2019), but also affects the flavor and shear force of pork (Zhang et al., 2021). Furthermore, consumers are more likely to accept meat with higher IMF (Malgwi et al., 2022). However, the molecular mechanism and transcriptional regulation process of IMF deposition are still unclear.

Alternative splicing (AS), a key step in the post-transcriptional regulation of gene expression (Lee and Rio, 2015), can directly participate in the pre-transcriptional and post-transcriptional processes and even the processing of RNA and affect the binding between downstream proteins and nucleic acids, proteins, or membranes, and change the localization and enzymatic properties of proteins (Kelemen et al., 2013). AS occurs in different stages of pig development and in different tissues and organs, and affects various economic traits of pig. AS is also involved in the regulation of biological processes and pathways, such as fatty acid degradation (Hao et al., 2022), glyceride metabolism (Li et al., 2012), glucose response (Liu et al., 2010), and fat deposition (Fang et al., 2021).

With the development of the third-generation sequencing (TGS) platform, the nanopore sequencing technology initiated by single-molecule real time (SMRT) technology and Oxford Nanopore Technologies (ONT) (Deamer et al., 2016) have shown great advantages in the identification of AS isoforms (Weirather et al., 2017; Beiki et al., 2019). However, an analysis of ONT-based AS for IMF development in pigs has not yet been conducted.

In this study, we used combined ONT and RNA sequencing analysis to identify potential differential AS factors affecting IMF deposition in pig longissimus dorsi muscle with high and low IMF content. In addition, key candidate genes and AS events affecting IMF deposition were screened by co-expression and PPI analysis. The results can facilitate the study of transcriptional regulatory processes and provide new insights into AS studies in livestock. At the same time, they can be used to provide a new method for improving IMF content in pig breeding.

Materials and Methods

Sample Collection and Determination of Intramuscular Fat Content

The animal materials used in this study were all from the F2 generation resource population of large white pig × Min pig at the Changping Experimental Base of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences. The pig raising and slaughtering environment, tissue collection, and IMF measuring methods were as described in our previous study (Wang et al., 2021). Thirty samples were screened for full-length transcript sequencing based on IMF content: fifteen samples with high IMF content and fifteen samples with low IMF content. The content in the high-IMF group ranged from 4.07% to 12.18%, and that in the low-IMF group ranged from 0.78% to 1.60%. The sample IMF content details are shown in Table 1 and Supplementary Table S13.

TABLE 1
www.frontiersin.org

TABLE 1. Sample IMF content details.

Oxford Nanopore Technologies and Strand-Specific RNA-Seq Library Construction and Sequencing

The ONT library construction was performed according to the standard protocol provided by the ONT platform. Before full-length transcript sequencing, 1 μg of total RNA was prepared for the cDNA library using a cDNA-PCR sequencing kit (SQK-PCS109). The libraries were then created using a sequencing library prep kit. We added the defined PCR adapters directly to both ends of the first-strand cDNA. The cDNA library was established by 20 cycles of PCR amplification with LongAmp Tag. The PCR products were then subjected to ONT adapter ligation using T4 DNA ligase. Agencourt XP beads were used for DNA purification according to the ONT protocol. The final cDNA library was added to the FLO-MIN109 flow cell, and the library was sequenced using the MinION Mk1B sequencer. Strand-specific RNA sequencing data were obtained to validate and quantify the sequencing data of the ONT full-length transcripts. Approximately 2 μg of RNA per sample was used to prepare Illumina sequencing libraries, and 30 prepared libraries were sequenced using the Illumina HiSeq 4,000 platform with a 150 bp pair-end sequencing strategy according to the manufacturer’s instructions.

Data Analysis of Nanopore Raw Data

For raw data from the Nanopore sequencing fast5 format, Guppy software in the MinKNOW2.2 software package was used to perform base calling and conversion to fastq format. Nanoplot (De Coster et al., 2018) was used to filter the converted data. The clean data were compared with the Silva rRNA database using sortmerna (Kopylova et al., 2012) to delete ribosomal RNA and for subsequent analysis. The sequence containing primers at both ends were determined to be a full-length sequence using pychopper (https://github.com/nanoporetech/pychopper/) software. The obtained full-length sequences were aligned with the pig reference genome (Sus scrofa 11.1.96) using minnimap2 software (Li, 2018), and the data were cleaned and clustered using pinfish (https://github.com/nanoporetech/pinfish/) to obtain consistency sequences. cDNA_Cupcake (https://github.com/Magdoll/cDNA_Cupcake/) was used to filter the sequences with identity ≤0.9 and coverage ≤0.85. Alignments with differences in only the 5′ end exons were merged, and 43,668 strips of gff file redundancies were obtained for the sequence information of the remaining transcripts. In addition, gffcompare (https://github.com/gpertea/gffcompare) was used to compare the transcripts obtained from full-length sequencing with the known transcripts of the genome, to discover new genes and transcripts and to supplement the genome annotation.

Identification of Alternative Splicing Events

To detect all AS events, we generated gff files after de-redundancy and used AStalavista software (Foissac and Sammeth, 2007) to perform AS analysis. Five types of AS events were obtained: exon skipping (ES), intron retention (IR), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), and mutually exclusive exons (MEE). The results of the AStalavista analysis were then transformed to obtain AS event information for all transcripts.

Validation of the Existence of Alternative Splicing Events by RNA Sanger Sequencing

To verify the accuracy of the identified alternative splicing events, we randomly selected several splicing isoforms for PCR validation. Specific primers upstream and downstream of the alternative splicing were designed using Primer 5. The PCR amplification reaction mixture was 25 μl, including 12.5 μl of Premix Taq (Takara RR902A), 1 μl of template cDNA, 1 μl of each primer, and 9.5 μl of sterile water. PCR products were examined by 1.0% agarose gel electrophoresis and sequenced using the Sanger sequencing method. The primers used in this study are listed in Supplementary Table S12.

Functional Annotation and Enrichment Analysis of Full-Length Transcripts

To obtain comprehensive functional information annotations for all transcripts, seven databases were searched using BLASTX v.2.2.26 (E-value < 1 × 10−5): Gene Ontology (GO, http://www.geneontology.org), Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/), Cluster of Orthologous Groups of proteins (COG, http://www.ncbi.nlm.nih.gov/COG/), Clusters of Orthologous Groups for eukaryotic complete genomes (KOG, ftp://ftp.ncbi.nih.gov/pub/COG/KOG/kyva), Non-Redundant Protein Sequences Database (NR, http://www.ncbi.nlm.nih.gov/), Swiss-Prot (http://www.expasy.org/sprot/), and Pfam database (Pfam, http://pfam.xfam.org/) (Finn et al., 2011) to annotate all de-redundant full-length transcripts. The R package clusterProfiler (Yu et al., 2012) was used to perform GO and KEGG pathway enrichment analysis on alternatively spliced genes, and the significant enrichment criterion was a q-value of <0.05.

Quantitative Gene Expression Analysis, Differentially Expressed Splicing Isoforms, and Differentially Expressed Genes Analysis

The full-length transcripts generated by ONT sequencing were annotated as reference genomes, and the expression levels of all transcripts were further calculated based on the data generated by Illumina sequencing. Hisat2 was used to estimate gene expression levels and to map clean data onto reference sequences. Stringtie with parameter -e was used to limit read alignment processing to estimate and output only assembled transcripts that matched a given reference transcript. We used python script prepDE.py to obtain read count files from the mapping results.

Only the expression data of alternatively spliced transcripts were extracted, and DESeq2 (Love et al., 2014) was used to analyze DESIs and DEGs between the high- and low-IMF groups. The DEseq2 package was used for analysis, the parameter design was set to types and conditions, the type parameter was the sequencing batch, and the conditions parameter was the test grouping. DESIs and DEGs were filtered according to the following criteria: padjust <0.05 and |log2 FoldChange| ≥ 1.

Expression Proportion Analysis and Protein Structure Prediction of Differentially Expressed Splicing Isoforms

In order to explore the regulatory effect of the identified DESIs on gene expression, we analyzed the proportion of all splicing isoforms in their gene expression:

Expression proportion of splicing isoforms (count%)=splicing isoform readsgene reads

Protein hydrophobicity (ProtScale, https://web.expasy.org/protscale/) and protein tertiary structure prediction (SWISS-MODEL, https://www.swissmodel.expasy.org/) of amino acid sequences of splicing isoforms were analyzed to explore the potential effect mechanism of AS on gene function.

Co-Expression and Protein-Protein Interaction Analysis

In this study, to explore the key factors of AS affecting IMF content, we also performed co-expression analysis using weighted gene co-expression network analysis (WGCNA) on all transcripts that were alternatively spliced. Splicing isoforms with a weighted threshold of 0.08 in the screening module were used as hub splicing isoforms. STRING v11.0 (https://string-db.org/) was used to analyze the protein interaction of splicing isoform genes in the module. Cytoscape (version 3.7.2) was used to display the protein interaction map, and the plugin cytoHubba based on degree values was used to identify hub genes in the PPI network.

Results

Overview of Oxford Nanopore Technologies and Strand-Specific RNA Sequencing Data

Full-length transcript sequencing was performed on 30 samples, and the clean data generated by each sample sequencing amounted to 5.52 GB on average. The average number of sequences per sample generated by nanopore sequencing was 7,290,257, and the quality values of the samples were all above Q8 (Table 2). A total of 148 GB of raw data was obtained by strand-specific RNA sequencing, and the annotation files obtained by ONT sequencing were aligned and assembled, with an average overall alignment rate of 95.78%. An overview of the next-generation sequencing data is shown in Supplementary Table S1.

TABLE 2
www.frontiersin.org

TABLE 2. Clean reads generated by Nanopore sequencing.

Overview of Full-Length Transcripts

For the 30 individuals, the average number of full-length sequences was 5,081,792 and the average matching proportion of full-length sequences was 79.67%, as shown in Supplementary Table S2. All consensus transcript sequences were aligned to the reference genome by minimap2 software and then subjected to de-redundancy analysis. Finally, 43,688 full-length transcripts were obtained. Mapping of the Sus scrofa genome by gffcompare identified gene loci and transcript isoforms. A total of 15,130 genes were detected; 10,808 (71.43%) known genes were matched and 4,322 (28.57%) genes were newly discovered (Figure 1A). A total of 43,688 transcripts were detected, among which only 6,086 (13.93%) transcripts were known. A total of 30,795 (70.49%) novel transcripts of known genes were identified, and 6,807 novel transcripts of novel genes were identified (15.58%) (Figure 1B). A total of 6,548 (43.27%) genes were identified with two or more transcripts (Figure 1C). Troponin t3, fast skeletal type (TNNT3) had the highest number of transcripts (143 novel and 11 known).

FIGURE 1
www.frontiersin.org

FIGURE 1. Overview of ONT transcripts. (A) Venn diagram of genes identified in ONT data and annotated genomes; (B) Distribution diagram of three types of transcripts in ONT data; (C) Distribution diagram of gene multi-transcripts; (D) Venn diagram of all new transcripts annotation result.

To elucidate their functions, all newly discovered transcripts were annotated to the GO, KEGG, COG, KOG, Pfam, NR, and SwissProt databases. A total of 29,539 transcripts could be annotated in at least one database (Figure 1D).

Alternative Splicing Identification

AStalavista was used to identify AS of de-redundant transcripts; a total of 10,949 splicing isoforms (25.06%, 43,688 transcripts detected) were generated, and 14,728 AS events were detected (Supplementary Table S4). The largest proportion of AS events was ES (5,040 events, 34.22%), and the smallest proportion was MEE (1,206 events, 8.19%; Figure 2A; Table 3). A total of 2,691 genes (17.79%, 15,130 genes detected by ONT) were alternatively spliced during transcription, and the top five genes with the most splicing isoforms were TNNT3 (139 splicing isoforms), troponin t1, slow skeletal type (TNNT1, 131 splicing isoforms), Y-box binding protein 3 (YBX3, 91 splicing isoforms), myosin binding protein C2 (MYBPC2, 72 splicing isoforms), and nebulin (NEB, 71 splicing isoforms) (Figure 2B).

FIGURE 2
www.frontiersin.org

FIGURE 2. Overview of alternative splicing. (A) The proportion of different AS events; (B) The number of different gene splicing isoforms; (C) Verify different AS events.

TABLE 3
www.frontiersin.org

TABLE 3. Number and percentage of AS events.

Validation of Selected Alternative Splicing Isoforms

To confirm the accuracy of the identified AS events, we randomly selected several and verified ONT.6273.5, ONT.15020.2, ONT.1066.3, ONT.8241.1, ONT.9799.3, and ONT.9654.4 for validation (Figure 2C). The results support the accuracy of the identified AS events.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Analysis of Genes With Specific Splicing Isoforms

After de-redundancy, we counted the number of de-redundant transcripts for each sample. De-redundant transcripts that were not present in the high- and low-IMF groups but were present in the other group were defined as low- or high-IMF group specific splicing isoforms. In total, 1,605 splicing isoforms were found only in the high-IMF group and 2,235 in the low-IMF group. GO and KEGG analysis was performed on genes of these differentially splicing isoforms using clusterProfiler (Supplementary Table S5). The results indicate that the significant biological process (BP) terms of high-IMF group genes were mainly ribonucleoprotein complex biogenesis, ribonucleoprotein complex assembly, and ribonucleoprotein complex subunit organization; the significant cellular component (CC) terms were ribonucleoprotein complex, proteasome complex, and endopeptidase complex; and the significant molecular function (MF) terms were RNA binding, structural constituent of ribosome, and structural molecule activity (Figure 3A). The enriched KEGG pathways were ribosome, spliceosome, and proteasome (Figure 3B). GO analysis of differentially splicing isoforms in the low-IMF group showed that significant the BP terms were mainly translation, peptide biosynthetic process, and cellular protein-containing complex assembly; significant CC terms were ribonucleoprotein complex, proteasome complex, and endopeptidase complex; and significant MF terms were RNA binding, structural constituent of ribosome, and structural molecule activity (Figure 3C). The enriched KEGG pathways were mainly ribosome, Parkinson disease, and Huntington disease (Figure 3D).

FIGURE 3
www.frontiersin.org

FIGURE 3. GO and KEGG analysis of splicing isoforms specifically occurring in H and L groups. (A) GO results of splicing isoforms that specifically occur in high-IMF group; (B) KEGG results of splicing isoforms that occur specifically in high-IMF group; (C) GO results of splicing isoforms that specifically occur in low-IMF group; (D) KEGG results of splicing isoforms that occur specifically in low-IMF group.

Differentially Expressed Splicing Isoforms, Differentially Expressed Genes, and Expression Proportion Analysis

Because the ONT method has a disadvantage in quantitating expression, we used RNA-seq data to adjust the expression of ONT sequencing data. Combining ONT with RNA-seq data, a total of 10,760 (98.27%) splicing isoform expressions were matched and assembled. Through differential analysis of the expression data of splicing isoforms occurring between high- and low-IMF groups, 98 DESIs (Figure 4A and Supplementary Figure S1), including 49 upregulated and 49 downregulated splicing isoforms, were investigated (Supplementary Table S6). Among all DESIs, SERPINE1 mRNA binding protein 1 (SERBP1) had four upregulated splicing isoforms (ONT.13291.2, ONT.13291.18, ONT.13291.21, ONT.13291.24); myosin light chain 1 (MYL1) had three upregulated splicing isoforms (ONT.5824.56, ONT.5824.110, ONT.5824.117); and TNNT3 had three downregulated splicing isoforms (ONT.7235.51, ONT.7235.143, ONT.7235.171). Interestingly, TNNT1 had two splicing isoforms, with ONT.12529.9 showing an upregulated trend and ONT.12529.241 a downregulated trend.

FIGURE 4
www.frontiersin.org

FIGURE 4. Volcano plot of DESIs and DEGs. (A) Volcano plot of DESIs, showing the top 20 splicing isoforms ID with the smallest q-values; (B) Volcano plot of DEGs, showing the top 15 genes ID with the smallest q-values.

Through differential analysis of gene expression data between the high- and low-IMF groups, 15 DEGs were identified, including 3 upregulated and 12 downregulated genes (Figure 4B; Supplementary Table S7). Combining the results of DEGs and DESIs, among the downregulated DEGs, ONT.12957 [interferon alpha inducible protein 6 (IFI6)] had one transcript (ONT.12957.1) identified in the DESIs results, and ONT.4294 [growth arrest and DNA damage inducible gamma (GADD45G)] had two transcripts (ONT 0.4294.3, ONT.4294.9) identified in the DESIs results.

As the ONT sequencing transcripts were long and could cover one gene, we used a criterion to screen the DESIs most likely to be effective, which was the proportion of isoform expression to the expression of all isoforms in the same gene above 10%. As shown in Table 4, we found that ONT.9523.3, ONT.3832.5, ONT.191.5, ONT.15153.3, ONT.11047.3, ONT.16074.8, ONT.10763.6, ONT.7120.12, ONT.7120.13, ONT.14524.8, and ONT.10349.9 had differences in the proportion of expression between high- and low-IMF pigs (p < 0.05).

TABLE 4
www.frontiersin.org

TABLE 4. Proportion of splicing isoforms expression in Group H and L.

Protein Structural Analysis of Splicing Isoforms

For most of the above 11 isoforms, their expression change could not result in significant expression difference from other transcription expression except ONT.15153.3 and ONT.10349.2, as shown in Supplementary Table S8. From the results, we found that increased expression of ONT.15153.3 will cause reduced expression of ONT.15153.1, and increased expression of ONT.10349.9 will cause reduced expression of ONT.10349.2. We then predicted the protein hydrophobicity and tertiary structure of the genes. We found that the hydrophobicity and tertiary structure of ADP-ribosyltransferase 1 (ART1) located in ONT.15153.3 had changed, as shown in Figures 5A,B. Thus, we inferred that the occurrence of AS in ART1 may have an impact on its protein function, thereby regulating the expression of its gene function.

FIGURE 5
www.frontiersin.org

FIGURE 5. Protein structure analysis of ONT.15153.1, ONT.15153.3, ONT.10349.2 and ONT.10349.9. (A) Protein hydrophobicity analysis and protein tertiary structure analysis of ONT.15153.1 and ONT.15153.3; (B) Protein hydrophobicity analysis and protein tertiary structure analysis of ONT.10349.2 and ONT.10349.9.

Co-Expression Studies of Splicing Isoforms

A co-expression module of the splicing isoforms and IMF deposition was constructed using the WGCNA package tool. By plotting a cluster plot of sample correlations, it could be seen that there were no obvious outlier samples (Figure 6A). We selected the value where the scale-free network graph structure R square reaches 0.95, that is, power = 6 (Figure 6B). By constructing a co-expression network, a total of 42 co-expression modules with different splicing isoforms were identified. Each module was displayed in a different color and contained different numbers of splicing isoforms. Hierarchical clustering and correlation heatmaps were constructed with the 42 modules and correlation heatmaps (Figures 6C,D). Correlation analysis was performed between the co-expression modules and IMF content traits, and correlations and p-values between different modules and traits were obtained (Figure 6E).

FIGURE 6
www.frontiersin.org

FIGURE 6. Co-expression analysis results. (A) Clustering dendrogram of 30 samples; (B) Determination of soft threshold power in co-expression analysis; (C) Hierarchical clustering dendrogram of all modules; (D) Correlation diagram between modules; (E) Relationship between all modules and IMF traits correlation graph.

Through co-expression analysis of splicing isoforms, we found a module (MEblack) that was significantly associated with IMF traits. We then plotted the significance of splicing isoform module membership for this module (Figure 7A). At the same time, we also performed correlation analysis on the splicing isoforms in the module and traits (Figure 7B).

FIGURE 7
www.frontiersin.org

FIGURE 7. Correlation analysis between MEblack and IMF. (A) Dot plot of correlation between splicing isoforms in MEblack and IMF; (B) Pheatmap of splicing isoforms in MEblack and 30 samples.

Module MEblack aggregated 191 splicing isoforms (Supplementary Table S9), among which 41 hub splicing isoforms were identified by setting a weighted threshold (Supplementary Table S10). Corresponding to 162 named genes, they were imported into the String database to construct a PPI network, and Cytoscape was used to display the PPI results (Figure 8A). Using the cytoHubba plugin, 30 hub genes were identified by the degree method (Figure 8B and Supplementary Table S11). The genes and splicing isoforms were mapped and compared with the hub splicing isoforms screened by the co-expression analysis results, and eight hub splicing isoforms were finally identified (ONT.2674.11, ONT.6603.3, ONT.7429.11, ONT.7429.12, ONT.10638.13, ONT.10638.14, ONT.13247.4, ONT.15079.18) (Figure 8C). The corresponding genes were found to be mitochondrial ribosomal protein L27 (MRPL27), AAR2 splicing factor (AAR2), glycogen phosphorylase, muscle associated (PYGM), proteasome 26S subunit ubiquitin receptor, non-ATPase 4 (PSMD4), sodium channel modifier 1 (SCNM1), and heterogeneous nuclear ribonucleoprotein D like (HNRNPDL).

FIGURE 8
www.frontiersin.org

FIGURE 8. Co-expression and PPI analysis results. (A) The gene-PPI network of splicing isoforms in MEblack; (B) The PPI network of the Degree top 30 gene; (C) The venn diagram of the hub splicing isoforms of MEblack and the hub splicing isoforms of PPI.

Discussion

Full-length transcript sequencing has been shown to be an effective method for studying the post-transcriptional regulation of genes (Cole et al., 2020), and full-length transcripts can be detected by end-to-end sequencing of most transcript sequences, thereby overcoming the errors caused by the complex transcript assembly steps (Ali et al., 2021) required for short-read sequencing. Using ONT sequencing, a total of 43,688 full-length transcripts were identified in this study, 4,322 novel genes and 30,795 novel transcripts. In longissimus dorsi samples, exon skipping and alternative 5′ splice site were the most frequent AS events, accounting for 34.22 and 28.75%, respectively. In a previous study, exon skipping and alternative 5′ splice site were the most common AS events in mammalian pre-mRNA (Sammeth et al., 2008), which is consistent with our study.

The analysis of AS events showed that the five genes with the most splicing isoforms were TNNT3, TNNT1, YBX3, MYBPC2, and NEB. TNNT3 and TNNT1 encode troponin T in fast and slow skeletal muscles, respectively, and play important roles in muscle contraction and regulation of calcium concentration (Wei and Jin, 2016). YBX3 plays an important role in binding to nucleic acids and transcription factors (Lyabin et al., 2020). MYBPC2 is involved in the regulation of actin binding, protein binding, muscle structural components, cell adhesion, and myofilament sliding (Chen et al., 2021). NEB is a structural component of sarcomeres that extend along filaments and regulate key functions of skeletal muscle (Gineste et al., 2020). In pig longissimus dorsi samples, genes related to skeletal muscle function produced numerous splicing isoforms, and skeletal muscle fibers were also found to be a factor affecting IMF deposition (Veloso et al., 2018), suggesting that AS plays an important role in muscle contraction and activity.

In our study, we also investigated 1,605 and 2,235 splicing isoforms generated only in samples from high-IMF and low-IMF groups, respectively. BP entries enriched in the high-IMF group were significantly different from BP entries enriched in the low-IMF group. BP entries enriched in the high-IMF group were ribonucleoprotein complex assembly, ribonucleoprotein complex subunit organization, and ribonucleoprotein complex biogenesis, while those in the low-IMF group were ubiquitin-dependent protein catabolic process, modification-dependent protein catabolic process, and modification-dependent macromolecule catabolic process. In the KEGG analysis, the significantly enriched pathways in the high- and low-IMF groups were all related to ribosomes, spliceosomes and diseases, glycolysis and gluconeogenesis, and fatty acid metabolism partly caused by abnormal muscle function. The functions of splicing isoforms in the high- and low-IMF groups are mostly related to the regulation of transcriptional processes and muscle function and a series of biological processes of fat synthesis and metabolism. The above results indicate that in the longissimus dorsi of pigs, AS is involved in the regulation of transcription and translation processes, and the spliceosome in the regulation of the occurrence of AS events. AS regulates protein synthesis and catabolism and other processes, which affect the differential functional expression of proteins. At the same time, AS also participates in the regulation of complex processes such as muscle function and fat anabolism.

The expression abundance of all full-length transcripts was obtained by comprehensive analysis of ONT and RNA-seq data. A total of 98 DESIs were identified in high- and low-IMF groups. Four upregulated DESIs were identified for SERBP1, which is highly expressed in adipose tissue and regarded as a signal of fat deposition (Duarte et al., 2016). SERBP1 produces plasminogen activator inhibitor 1 RNA-binding protein (PAIRBP1), an inflammatory adipokine mainly present in adipocytes. Studies have shown that fat deposition can lead to steatosis, triggering the development of fibrosis and steatohepatitis, resulting in upregulation of PAIRBP1 expression (Ota et al., 2007). In addition, the expression of SERBP1 can inhibit angiostatin to promote angiogenesis (Nishioka et al., 2011), and the processes of adipogenesis and angiogenesis are synchronized. Our study identified four DESIs of SERBP1 in the high-IMF group, leading us to speculate that AS may play an important role in IMF deposition by regulating SERBP1 expression. Moreover, three upregulated DESIs were identified for MYL1, which encodes myosin light chain 1. Studies have shown that MYL3, the paralogous gene of MYL1, plays a negative role in adipogenic differentiation (Zhang et al., 2010). MYL1 is inversely correlated with muscle tenderness (Gagaoua et al., 2020) and is highly expressed in tougher meats and in pigs with leaner carcasses. MYL1 is present in fast-twitch fibers, which contain less myoglobin, and ATP is produced by glycolytic metabolism, resulting in more fast-twitch fibers and lower IMF content in tougher meat (Rodrigues et al., 2017). In our study, multiple upregulated DESIs of MYL1 were identified in samples from the high-IMF group, and we speculated that AS may produce functionally opposite splicing isoforms of MYL1 to make it function abnormally, and these isoforms may compete for MYL1 translation substrates and reduce MYL1 expression in high-IMF samples. This result suggests that AS may enhance gene function or generate functionally opposite splicing isoforms to alter gene expression and function.

Three downregulated DESIs were identified for TNNT3, which is involved in encoding fast skeletal muscle fibers, and whose dysfunction contributes to the development of muscular dystrophy symptoms. Castrated male bulls produce more IMF, and TNNT3 showed a trend of upregulation after castration (Silva et al., 2019). The promotion of IMF deposition by TNNT3 has been verified in chicken (Liu et al., 2016) and cattle (Bong et al., 2009). Interestingly, a study of the DNA-binding motif of TNNT3 found that its expression levels were inversely correlated with whole-body fat mass (rTNNT3 ∼ fat mass = −0.49) (Nunez Lopez et al., 2018). This suggests that different splicing isoforms of TNNT3 play different roles in the process of IMF deposition, which was also confirmed (Schilder et al., 2011). We found in DESIs analysis that TNNT1 had one upregulated and one downregulated splicing isoform. TNNT1 is involved in encoding slow skeletal muscle fibers. Studies have found that TNNT1 is directly associated with obesity traits (Pierzchala et al., 2014), and increased TNNT1 expression is positively correlated with triglycerides (Nayak et al., 2010). TNNT1 showed a downward trend in castrated bulls (Silva et al., 2019). The different functions of TNNT1 may be caused by the opposite functions of different splicing isoforms, which is consistent with the results of our analysis. From the multiple regulatory roles of AS of MYL1, TNNT3, and TNNT1, it can be seen that the splicing patterns of skeletal muscle are diverse and complex, and AS has a critical impact on the development and function of muscle fibers and fat deposition. The expression abundance of different splicing isoforms may promote or reverse the realization of gene function, which greatly enriches proteome diversity.

Combining the results of DEGs and DESIs, it was found that the IFI6 gene and its splicing isoform ONT.12957.1 were both downregulated in the low-IMF group. The GADD45G gene and its splicing isoforms ONT.4294.3 and ONT.4294.9 were also downregulated in the low-IMF group. IFI6 is a gene that may play an important role in meat quality and carcass indicators and is regarded as a candidate gene for improved pork quality (Kayan et al., 2011). GADD45G is an inflammation-related gene associated with oxidative stress and cytokine secretion signaling pathways. The expression of GADD45G is upregulated in the liver in high-fat-diet-induced non-alcoholic fatty liver disease (NAFLD), and its expression decreases when the high-fat condition is improved (Wang et al., 2022), which may be related to its effect of on hepatic fibrosis and the occurrence of chemical transformation (Hong et al., 2016). The above results suggest that differential expression of spliced isoforms may cause differences in gene expression.

By analyzing the proportion of expression of each DESI in all transcripts, we explored the important role of splicing isoforms in the realization of gene function. Among the DESIs in the high- and low-IMF groups, we found that the average expression proportion of ONT.9523.3, ONT.3832.5, ONT.15153.3, ONT.11047.3, ONT.16074.8, ONT.10763.6, ONT.7120.12, ONT.7120.13, ONT.14524.8, ONT.10349.9, and ONT.191.5 in samples from the two groups was >10%, and the difference between the two was significant (p < 0.05). Then, we compared the expression ratio of splicing isoforms and other transcripts of the gene. The proportion of ONT.15153.1 expression in the low-IMF group was significantly lower than that in the high-IMF group. We speculated that increased expression of ONT.15153.3 would cause decreased expression of ONT.15153.1. At the same time, increased expression of ONT.10349.9 also caused decreased expression of ONT.10349.2. The analysis of its protein hydrophobicity and tertiary structure showed that the occurrence of AS events will lead to different coding information of transcripts, which will change the hydrophobic region of the amino acid sequence, and then affect the conformation of the tertiary structure of the protein. The gene encoding ONT.15153.3 is ART1, which catalyzes the ADP-ribosylation of arginine residues in the protein. It is highly expressed in skeletal muscle and is associated with the formation of myotubes and muscle fibers (Leutert et al., 2018). The gene encoding ONT.10349.9 is RAB2A, a member of the RAS oncogene family (RAB2A) and a membrane-bound protein involved in vesicle fusion and trafficking. RAB2A knockdown inhibits glucose-stimulated insulin secretion (Sugawara et al., 2014). According to the above results, the occurrence of AS may change the expression of some transcripts of the same gene, thereby affecting gene function.

Through differential analysis of the expression abundance of alternatively spliced transcripts, it was found that AS plays a variety of regulatory roles in the process of fat deposition. In addition, by constructing a co-expression network, the splicing isoforms of eight centers common to the co-expression and PPI networks were screened, which belonged to six genes: MRPL27, AAR2, PYGM, PSMD4, SCNM1, and HNRNPDL. MRPL27 encodes a mitochondrial ribosomal protein, and its abnormal function may cause small steric dysfunction and muscle atrophy (Bogatikov et al., 2020). MRPL27 is also regarded as an NAFLD-related hub gene (Zeng et al., 2021). AAR2 is a splicing factor that is a component of the U5 snRNP complex and is required for spliceosome assembly and pre-mRNA splicing (Weber et al., 2013). PYGM encodes glycogen phosphorylase, which provides energy for muscle contraction, and is associated with glycogen storage in skeletal muscle (Tarnopolsky, 2018). PSMD4 encodes proteasome 26S subunit ubiquitin receptor, which is a proteasome-related gene that is significantly downregulated in hypercholesterolemia (Loke et al., 2017). SCNM1 encodes a zinc finger protein and putative splicing factor, and its mis-splicing will lead to serious lesions of sodium channels in mice (Buchner et al., 2003). HNRNPDL encodes heterogenetic nuclear ribonucleoprotein d like, which is an important transcription regulator that regulates the AS of hundreds of genes (Li et al., 2019).

Conclusion

Our study shows that 17.79% of the genes in the longissimus dorsi muscle produced splicing isoforms. Specific splicing isoforms in the high- and low-IMF groups were found to be related to the regulation of the transcription process and muscle fiber function. MYL1, TNNT3, and TNNT1 produced different splicing isoforms through AS, which can promote or slow down the process of IMF deposition. The splicing isoforms of IFI6 and GADD45G were differentially expressed, which may cause differences in gene expression, and the expression of the splicing isoforms of ART1 and RAB2A caused changes in the expression of other transcripts. Several genes of hub splicing isoforms, such as MRPL27, AAR2, PYGM, PSMD4, SCNM1, and HNRNPDL, were also investigated with regard to iIMF. These genes may be associated with IMF deposition, but there is currently no direct evidence and further studies are needed to explore their function. Our study preliminarily shows that AS may affect fat deposition in many forms, which also provides further insight into the mechanism of AS regulating IMF content. Our results identified many candidate splice isoforms in IMF that can be used as molecular markers in the breeding and improvement of pig IMF content after validation.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The animal study was reviewed and approved by Animal Ethics Committee of the Institute of Animal Science, Chinese Academy of Agricultural Sciences.

Author Contributions

LgW and LxW conceived and designed the experimental plan, and ZS participated in the bioinformatics analysis and drafted the manuscript. LgW, JW, XH, HY and LxW helped revise this manuscript. All authors read and approved the final manuscript.

Funding

This work was jointly supported by the National Natural Science Foundation of China (Grant Nos 31872337 and 31501919), China Agriculture Research System of MOF and MARA (CARS-35), and Agricultural Science and Technology Innovation Project (ASTIP-IAS02). The funding bodies had no role in the design of the study, or collection, analysis and interpretation of data and writing the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank the researchers at our laboratories for their dedication and hard work.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.877646/full#supplementary-material

Supplementary Figure S1 | Cluster analysis of all samples with differentially splicing isoforms.

References

Alfaia, C. M., Lopes, P. A., Madeira, M. S., Pestana, J. M., Coelho, D., Toldrá, F., et al. (2019). Current Feeding Strategies to Improve Pork Intramuscular Fat Content and its Nutritional Quality. Adv. Food Nutr. Res. 89, 53–94. doi:10.1016/bs.afnr.2019.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Ali, A., Thorgaard, G. H., and Salem, M. (2021). PacBio Iso-Seq Improves the Rainbow Trout Genome Annotation and Identifies Alternative Splicing Associated with Economically Important Phenotypes. Front. Genet. 12, 683408. doi:10.3389/fgene.2021.683408

PubMed Abstract | CrossRef Full Text | Google Scholar

Beiki, H., Liu, H., Huang, J., Manchanda, N., Nonneman, D., Smith, T. P. L., et al. (2019). Improved Annotation of the Domestic Pig Genome through Integration of Iso-Seq and RNA-Seq Data. BMC Genomics 20 (1), 344. doi:10.1186/s12864-019-5709-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Bogatikov, E., Lindblad, I., Punga, T., and Punga, A. R. (2020). miR-1933-3p Is Upregulated in Skeletal Muscles of MuSK+ EAMG Mice and Affects Impa1 and Mrpl27. Neurosci. Res. 151, 46–52. doi:10.1016/j.neures.2019.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Bong, J., Cho, K., and Baik, M. (2009). Comparison of Gene Expression Profiling between Bovine Subcutaneous and Intramusular Adipose Tissues by Serial Analysis of Gene Expression. Cell Biol. Int. 34 (1), 125–133. doi:10.1042/CBI20090046

PubMed Abstract | CrossRef Full Text | Google Scholar

Buchner, D. A., Trudeau, M., and Meisler, M. H. (2003). SCNM1, a Putative RNA Splicing Factor that Modifies Disease Severity in Mice. Science 301 (5635), 967–969. doi:10.1126/science.1086187

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Z., Li, X.-y., Guo, P., and Wang, D.-l. (2021). MYBPC2 and MYL1 as Significant Gene Markers for Rhabdomyosarcoma. Technol. Cancer Res. Treat. 20, 153303382097966. doi:10.1177/1533033820979669

CrossRef Full Text | Google Scholar

Cole, C., Byrne, A., Adams, M., Volden, R., and Vollmers, C. (2020). Complete Characterization of the Human Immune Cell Transcriptome Using Accurate Full-Length cDNA Sequencing. Genome Res. 30 (4), 589–601. doi:10.1101/gr.257188.119

PubMed Abstract | CrossRef Full Text | Google Scholar

De Coster, W., D’Hert, S., Schultz, D. T., Cruts, M., and Van Broeckhoven, C. (2018). NanoPack: Visualizing and Processing Long-Read Sequencing Data. Bioinformatics 34 (15), 2666–2669. doi:10.1093/bioinformatics/bty149

PubMed Abstract | CrossRef Full Text | Google Scholar

Deamer, D., Akeson, M., and Branton, D. (2016). Three Decades of Nanopore Sequencing. Nat. Biotechnol. 34 (5), 518–524. doi:10.1038/nbt.3423

PubMed Abstract | CrossRef Full Text | Google Scholar

Duarte, J. L. G., Cantet, R. J. C., Rubio, Y. L. B., Bates, R. O., Ernst, C. W., Raney, N. E., et al. (2016). Refining Genomewide Association for Growth and Fat Deposition Traits in an F2 Pig Population1. J. Anim. Sci. 94 (4), 1387–1397. doi:10.2527/jas.2015-0182

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, X., Xia, L., Yu, H., He, W., Bai, Z., Qin, L., et al. (2021). Comparative Genome-wide Alternative Splicing Analysis of Longissimus Dorsi Muscles between Japanese Black (Wagyu) and Chinese Red Steppes Cattle. Front. Vet. Sci. 8, 634577. doi:10.3389/fvets.2021.634577

PubMed Abstract | CrossRef Full Text | Google Scholar

Finn, R. D., Clements, J., and Eddy, S. R. (2011). HMMER Web Server: Interactive Sequence Similarity Searching. Nucleic Acids Res. 39, W29–W37. Web Server issue. doi:10.1093/nar/gkr367

PubMed Abstract | CrossRef Full Text | Google Scholar

Foissac, S., and Sammeth, M. (2007). ASTALAVISTA: Dynamic and Flexible Analysis of Alternative Splicing Events in Custom Gene Datasets. Nucleic Acids Res. 35, W297–W299. Web Server issue. doi:10.1093/nar/gkm311

PubMed Abstract | CrossRef Full Text | Google Scholar

Gagaoua, M., Bonnet, M., and Picard, B. (2020). Protein Array-Based Approach to Evaluate Biomarkers of Beef Tenderness and Marbling in Cows: Understanding of the Underlying Mechanisms and Prediction. Foods 9 (9), 1180. doi:10.3390/foods9091180

PubMed Abstract | CrossRef Full Text | Google Scholar

Gineste, C., Ogier, A. C., Varlet, I., Hourani, Z., Bernard, M., Granzier, H., et al. (2020). In Vivo characterization of Skeletal Muscle Function in Nebulin‐deficient Mice. Muscle Nerve 61 (3), 416–424. doi:10.1002/mus.26798

PubMed Abstract | CrossRef Full Text | Google Scholar

Hao, W., Yang, Z., Sun, Y., Li, J., Zhang, D., Liu, D., et al. (2022). Characterization of Alternative Splicing Events in Porcine Skeletal Muscles with Different Intramuscular Fat Contents. Biomolecules 12 (2), 154. doi:10.3390/biom12020154

PubMed Abstract | CrossRef Full Text | Google Scholar

Hong, L., Sun, Q.-F., Xu, T.-Y., Wu, Y.-H., Zhang, H., Fu, R.-Q., et al. (2016). New Role and Molecular Mechanism of Gadd45a in Hepatic FIbrosis. Wjg 22 (9), 2779–2788. doi:10.3748/wjg.v22.i9.2779

PubMed Abstract | CrossRef Full Text | Google Scholar

Kayan, A., Uddin, M. J., Cinar, M. U., Grosse-Brinkhaus, C., Phatsara, C., Wimmers, K., et al. (2011). Investigation on Interferon Alpha-Inducible Protein 6 (IFI6) Gene as a Candidate for Meat and Carcass Quality in Pig. Meat Sci. 88 (4), 755–760. doi:10.1016/j.meatsci.2011.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Kelemen, O., Convertini, P., Zhang, Z., Wen, Y., Shen, M., Falaleeva, M., et al. (2013). Function of Alternative Splicing. Gene 514 (1), 1–30. doi:10.1016/j.gene.2012.07.083

PubMed Abstract | CrossRef Full Text | Google Scholar

Kopylova, E., Noé, L., and Touzet, H. (2012). SortMeRNA: Fast and Accurate Filtering of Ribosomal RNAs in Metatranscriptomic Data. Bioinformatics 28 (24), 3211–3217. doi:10.1093/bioinformatics/bts611

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, Y., and Rio, D. C. (2015). Mechanisms and Regulation of Alternative Pre-mRNA Splicing. Annu. Rev. Biochem. 84, 291–323. doi:10.1146/annurev-biochem-060614-034316

PubMed Abstract | CrossRef Full Text | Google Scholar

Leutert, M., Menzel, S., Braren, R., Rissiek, B., Hopp, A.-K., Nowak, K., et al. (2018). Proteomic Characterization of the Heart and Skeletal Muscle Reveals Widespread Arginine ADP-Ribosylation by the ARTC1 Ectoenzyme. Cell Rep. 24 (7), 1916–1929. doi:10.1016/j.celrep.2018.07.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H. (2018). Minimap2: Pairwise Alignment for Nucleotide Sequences. Bioinformatics 34 (18), 3094–3100. doi:10.1093/bioinformatics/bty191

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R. Z., Hou, J., Wei, Y., Luo, X., Ye, Y., and Zhang, Y. (2019). hnRNPDL Extensively Regulates Transcription and Alternative Splicing. Gene 687, 125–134. doi:10.1016/j.gene.2018.11.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Suh, Y., Kim, E., Moeller, S. J., and Lee, K. (2012). Alternative Splicing and Developmental and Hormonal Regulation of Porcine Comparative Gene Identification-58 (CGI-58) mRNA1. J. Anim. Sci. 90 (12), 4346–4354. doi:10.2527/jas.2012-5151

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Fu, R., Liu, R., Zhao, G., Zheng, M., Cui, H., et al. (2016). Protein Profiles for Muscle Development and Intramuscular Fat Accumulation at Different Post-Hatching Ages in Chickens. PLoS One 11 (8), e0159722. doi:10.1371/journal.pone.0159722

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, K., Wang, G., Zhao, S. H., Liu, B., Huang, J. N., Bai, X., et al. (2010). Molecular Characterization, Chromosomal Location, Alternative Splicing and Polymorphism of Porcine GFAT1 Gene. Mol. Biol. Rep. 37 (6), 2711–2717. doi:10.1007/s11033-009-9805-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Loke, S.-Y., Wong, P. T.-H., and Ong, W.-Y. (2017). Global Gene Expression Changes in the Prefrontal Cortex of Rabbits with Hypercholesterolemia And/or Hypertension. Neurochem. Int. 102, 33–56. doi:10.1016/j.neuint.2016.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2. Genome Biol. 15 (12), 550. doi:10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Lyabin, D. N., Eliseeva, I. A., Smolin, E. A., Doronin, A. N., Budkina, K. S., Kulakovskiy, I. V., et al. (2020). YB-3 Substitutes YB-1 in Global mRNA Binding. RNA Biol. 17 (4), 487–499. doi:10.1080/15476286.2019.1710050

PubMed Abstract | CrossRef Full Text | Google Scholar

Malgwi, I. H., Halas, V., Grünvald, P., Schiavon, S., and Jócsák, I. (2022). Genes Related to Fat Metabolism in Pigs and Intramuscular Fat Content of Pork: A Focus on Nutrigenetics and Nutrigenomics. Animals 12 (2), 150. doi:10.3390/ani12020150

PubMed Abstract | CrossRef Full Text | Google Scholar

Nayak, S. B., Pinto Pereira, L. M., Boodoo, S., Kimberlyali, A., Baptiste, C., Maraj, S., et al. (2010). Association of Troponin T and Altered Lipid Profile in Patients Admitted with Acute Myocardial Infarction. Arch. Physiol. Biochem. 116 (1), 21–27. doi:10.3109/13813450903397638

PubMed Abstract | CrossRef Full Text | Google Scholar

Nishioka, N., Matsuoka, T., Yashiro, M., Hirakawa, K., Olden, K., and Roberts, J. D. (2011). Linoleic Acid Enhances Angiogenesis through Suppression of Angiostatin Induced by Plasminogen Activator Inhibitor 1. Br. J. Cancer 105 (11), 1750–1758. doi:10.1038/bjc.2011.434

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunez Lopez, Y. O., Messi, M. L., Pratley, R. E., Zhang, T., and Delbono, O. (2018). Troponin T3 Associates with DNA Consensus Sequence that Overlaps with P53 Binding Motifs. Exp. Gerontol. 108, 35–40. doi:10.1016/j.exger.2018.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Ota, T., Takamura, T., Kurita, S., Matsuzawa, N., Kita, Y., Uno, M., et al. (2007). Insulin Resistance Accelerates a Dietary Rat Model of Nonalcoholic Steatohepatitis. Gastroenterology 132 (1), 282–293. doi:10.1053/j.gastro.2006.10.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Pierzchala, M., Hoekman, A. J. W., Urbanski, P., Kruijt, L., Kristensen, L., Young, J. F., et al. (2014). Validation of Biomarkers for Loin Meat Quality (M. Longissimus ) of Pigs. J. Anim. Breed. Genet. 131 (4), 258–270. doi:10.1111/jbg.12081

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodrigues, R. T. d. S., Chizzotti, M. L., Vital, C. E., Baracat-Pereira, M. C., Barros, E., Busato, K. C., et al. (2017). Differences in Beef Quality between Angus (Bos taurus Taurus) and Nellore (Bos taurus Indicus) Cattle through a Proteomic and Phosphoproteomic Approach. PLoS One 12 (1), e0170294. doi:10.1371/journal.pone.0170294

PubMed Abstract | CrossRef Full Text | Google Scholar

Sammeth, M., Foissac, S., and Guigó, R. (2008). A General Definition and Nomenclature for Alternative Splicing Events. Plos Comput. Biol. 4 (8), e1000147. doi:10.1371/journal.pcbi.1000147

PubMed Abstract | CrossRef Full Text | Google Scholar

Schilder, R. J., Kimball, S. R., Marden, J. H., and Jefferson, L. S. (2011). Body Weight-dependent Troponin T Alternative Splicing Is Evolutionarily Conserved from Insects to Mammals and Is Partially Impaired in Skeletal Muscle of Obese Rats. J. Exp. Biol. 214 (9), 1523–1532. doi:10.1242/jeb.051763

PubMed Abstract | CrossRef Full Text | Google Scholar

Silva, L. H. P., Rodrigues, R. T. S., Assis, D. E. F., Benedeti, P. D. B., Duarte, M. S., and Chizzotti, M. L. (2019). Explaining Meat Quality of Bulls and Steers by Differential Proteome and Phosphoproteome Analysis of Skeletal Muscle. J. Proteomics 199, 51–66. doi:10.1016/j.jprot.2019.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugawara, T., Kano, F., and Murata, M. (2014). Rab2A Is a Pivotal Switch Protein that Promotes Either Secretion or ER-Associated Degradation of (Pro)insulin in Insulin-Secreting Cells. Sci. Rep. 4, 6952. doi:10.1038/srep06952

PubMed Abstract | CrossRef Full Text | Google Scholar

Tarnopolsky, M. A. (2018). Myopathies Related to Glycogen Metabolism Disorders. Neurotherapeutics 15 (4), 915–927. doi:10.1007/s13311-018-00684-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Veloso, R., Lopes, P., Duarte, M., Silva, F., Saraiva, A., Chizzotti, M., et al. (2018). Research Article Expression of Lipid Metabolism and Myosin Heavy Chain Genes in Pigs Is Affected by Genotype and Dietary Lysine. Genet. Mol. Res. 17 (2). doi:10.4238/gmr16039904

CrossRef Full Text | Google Scholar

Wang, L., Zhou, Z.-Y., Zhang, T., Zhang, L., Hou, X., Yan, H., et al. (2021). IRLnc: a Novel Functional Noncoding RNA Contributes to Intramuscular Fat Deposition. BMC Genomics 22 (1), 95. doi:10.1186/s12864-020-07349-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, R., Yao, L., Lin, X., Hu, X., and Wang, L. (2022). Exploring the Potential Mechanism of Rhodomyrtus Tomentosa (Ait.) Hassk Fruit Phenolic Rich Extract on Ameliorating Nonalcoholic Fatty Liver Disease by Integration of Transcriptomics and Metabolomics Profiling. Food Res. Int. 151, 110824. doi:10.1016/j.foodres.2021.110824

PubMed Abstract | CrossRef Full Text | Google Scholar

Weber, G., Cristão, V. F., Santos, K. F., Jovin, S. M., Heroven, A. C., Holton, N., et al. (2013). Structural Basis for Dual Roles of Aar2p in U5 snRNP Assembly. Genes Dev. 27 (5), 525–540. doi:10.1101/gad.213207.113

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, B., and Jin, J.-P. (2016). TNNT1, TNNT2, and TNNT3: Isoform Genes, Regulation, and Structure-Function Relationships. Gene 582 (1), 1–13. doi:10.1016/j.gene.2016.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Weirather, J. L., de Cesare, M., Wang, Y., Piazza, P., Sebastiano, V., Wang, X.-J., et al. (2017). Comprehensive Comparison of Pacific Biosciences and Oxford Nanopore Technologies and Their Applications to Transcriptome Analysis. F1000Res 6, 100. doi:10.12688/f1000research.10571.210.12688/f1000research.10571.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16 (5), 284–287. doi:10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, F., Shi, M., Xiao, H., and Chi, X. (2021). WGCNA-based Identification of Hub Genes and Key Pathways Involved in Nonalcoholic Fatty Liver Disease. Biomed. Res. Int. 2021, 1–16. doi:10.1155/2021/5633211

CrossRef Full Text | Google Scholar

Zhang, L., Li, F., Guo, Q., Duan, Y., Wang, W., Yang, Y., et al. (2021). Balanced Branched‐chain Amino Acids Modulate Meat Quality by Adjusting Muscle Fiber Type Conversion and Intramuscular Fat Deposition in Finishing Pigs. J. Sci. Food Agric. doi:10.1002/jsfa.11728

CrossRef Full Text | Google Scholar

Zhang, Q., Lee, H.-G., Han, J.-A., Kim, E. B., Kang, S. K., Yin, J., et al. (2010). Differentially Expressed Proteins during Fat Accumulation in Bovine Skeletal Muscle. Meat Sci. 86 (3), 814–820. doi:10.1016/j.meatsci.2010.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Glossary

AS alternative splicing

IMF intramuscular fat

PPI protein-protein interaction

ONT Oxford Nanopore Technologies

PacBio Pacific Biosciences

TGS third-generation sequencing

SMRT single-molecule real time

WGCNA weighted gene co-expression network analysis

PCR polymerase chain reaction

ES exon skipping

IR intron retention

A5SS alternative 5′ splice site

A3SS alternative 3′ splice site

MEE mutually exclusive exons

GO Gene Ontology

KEGG Kyoto Encyclopedia of Genes and Genomes

COG Cluster of Orthologous Groups of proteins

KOG Clusters of orthologous groups for eukaryotic complete genomes

NR NCBI Non-redundant Protein Sequences Database

Pfam Pfamdatabase

BP biological process

CC cellular component

MF molecular function

DEGs differentially expressed genes

DESIs differentially expressed splicing isoforms

TNNT3 troponin t3, fast skeletal type

TNNT1 troponin t1, slow skeletal type

YBX3 Y-box binding protein 3

MYBPC2 myosin binding protein C2

NEB nebulin

SERBP1 SERPINE1 mRNA binding protein 1

MYL1 myosin light chain 1

IFI6 interferon alpha inducible protein 6

GADD45G growth arrest and DNA damage inducible gamma

ART1 ADP-ribosyltransferase 1

MRPL27 mitochondrial ribosomal protein L27

AAR2 AAR2 splicing factor

PYGM glycogen phosphorylase, muscle associated

PSMD4 proteasome 26S subunit ubiquitin receptor, non-ATPase 4

SCNM1 sodium channel modifier 1

HNRNPDL heterogeneous nuclear ribonucleoprotein D like

PAIRBP1 plasminogen activator inhibitor 1 RNA-binding protein

NAFLD nonalcoholic fatty liver disease

Keywords: alternative splicing, pig, intramuscular fat, nanopore sequencing, RNA-seq

Citation: Shu Z, Wang L, Wang J, Zhang L, Hou X, Yan H and Wang L (2022) Integrative Analysis of Nanopore and Illumina Sequencing Reveals Alternative Splicing Complexity in Pig Longissimus Dorsi Muscle. Front. Genet. 13:877646. doi: 10.3389/fgene.2022.877646

Received: 17 February 2022; Accepted: 22 March 2022;
Published: 11 April 2022.

Edited by:

Xiao Li, Northwest A&F University, China

Reviewed by:

Haiming Ma, Hunan Agricultural University, China
Paolo Zambonelli, University of Bologna, Italy

Copyright © 2022 Shu, Wang, Wang, Zhang, Hou, Yan and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Ligang Wang, ligwang@126.com

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.