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.

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.

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.

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.

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.

To obtain comprehensive functional information annotations for all transcripts, seven databases were searched using BLASTX v.2.2.26 (E-value < 1 × 10⁻⁵): 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.

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 |log₂ FoldChange| ≥ 1.

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 reads gene 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.

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.

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.

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).

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).

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).

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.

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).

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.

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).

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.

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).

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).

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).