Animal performance and longissimus lumborum meat quality and composition
These results have been reported in a companion paper34 and are mentioned here for contextualization purposes. Growth performance, carcass characteristics and meat quality traits were not significantly affected by dietary treatments. Pigs fed with C. vulgaris diets significantly increased total carotenoid deposition. For instance, CV + R increased the total chlorophylls and carotenoids by 2.1-fold compared to control. Microalga diets significantly increased n-3 PUFA and reduced the n-6:n-3 ratio in meat, thanks to increases in n-3 PUFA including 18:3n-3, 20:3n-3, 20:5n-3, 22:5n-3, and 22:6n-3.
Sequencing output and identification of expressed transcripts in the pig muscle transcriptome
In this study, we performed the analysis of the whole muscle transcriptome of each finishing pig fed with the following experimental diets: Control, CV, CV + R and CV + M. The total number of reads generated by RNA-Seq are reported in Supplementary Table S1. The average number of raw reads per experimental group was 65,140,430 for control, 38,830,084 for CV, 31,702,033 for CV + R and 30,641,589 for CV + M. Raw reads were trimmed, generating the average numbers of reads per sample for each experimental group of 45,219,804.5 (69.4%), 34,349,352.17 (88.5%), 30,768,719.67 (97.1%) and 27,093,411 (88.4%) for Control, CV, CV + R and CV + M, respectively. The trimmed reads were used for further analysis.
The trimmed reads were mapped in the Sus scrofa genome using hisat2 tool. The total average number of mapped reads was 40,839,191 (89.5%) for control, 30,977,391 (90%) for CV, 28,155,015 (91.5%) for CV + R and 24,558,272 (90.7%) for CV + M, while the average of uniquely mapped reads was 34,293,067 (74.7%) for control, 28,293,700 (81.3%) for CV, 25,692,735 (83.1%) for CV + R, and 21,706,256 (79.7%) for CV + M. The density of the mapped reads on different regions of the genome is displayed in Supplementary Figure S1.
The gene expression was quantified using the featureCounts tool, which reports the raw counts of reads that map to a single location in each of the exons (feature) that belong to that gene35. The top ten most abundantly expressed coding genes of muscle transcriptome in all experimental groups (raw counts from 340,811 to 1,219,730), ranked by absolute abundance were: ATP synthase F0 subunit 6 (ATP6), Cytochrome c oxidase subunit III (COX3), Nebulin (NEB), Enolase 3 (ENO3), NADH dehydrogenase subunit 4 (ND4), Myosin heavy chain 1 skeletal muscle (MYH1), Lactate dehydrogenase A (LDHA), Actin alpha 1 skeletal muscle (ACTA1), Myosin light chain phosphorylable fast skeletal muscle (MYLPF), and Myosin light chain 1 (MYL1) (Supplementary Fig. S2).
Transcriptome profile and differential expression
Among the muscle transcriptome data, we identified 781, 596, 744 and 575 genes, including protein-coding genes and ncRNA genes, that were expressed in the Control, CV, CV + R and CV + M groups, respectively. We also identified 17,811 genes that were expressed in all experimental groups and 24,584 genes that were expressed, at least, in one experimental group.
The DESeq2 tool was used to test for differential gene expression analysis between the following comparisons among experimental diets: Control vs CV, Control vs CV + R, and Control vs CV + M. In the pool of annotated genes from pig muscle transcriptome, we identified 1190 significant DEGs (245 upregulated and 945 downregulated) in the Control vs CV, 969 significant DEGs (277 upregulated and 692 downregulated) in the Control vs CV + R, and 3 significant DEGs (2 upregulated and 1 downregulated) in the Control vs CV + M. The number of non-redundant DEGs is 649 for Control vs CV and 428 for Control vs CV + R, with 538 DEGs common to these two comparisons and 3 DEGs common to the three comparisons. This is represented in the Venn diagram (Fig. 1). Based on the results from the DEGs analysis, a Multi-dimensional scaling (MDS) plot was performed in order to assess how the samples are distributed in a two-dimensional analysis according to leading logFC, and also a hierarchical clustering of treatment condition samples based on Pearson and Spearman correlations. A volcano plot of DEGs, which shows the relationship between FC and evidence of differential expression (-log FDR) was also developed.
Venn diagram with the distribution of the DEGs in the three experimental diet comparisons.
In the comparison Control vs CV, we observed in the MDS plot a good discrimination between the two experimental groups with samples from Control to occupy higher values of leading logFC dimension 1 with the exception of the sample 36S (Fig. 2A), which is according with the correlations (Fig. 2C). The red dots in the volcano plot highlight some of the identified DEGs (Fig. 2B).
(A) Multi-dimensional scaling (MDS) plot showing the relation between Control and CV diet samples. (B) Volcano plot showing the relationship between FC and evidence of differential expression (-log FDR) for Control vs CV. (C) Clustering of sample to sample correlations—Control vs CV. Hierarchical clustering of treatment condition samples based on Pearson (left) and Spearman correlations (right).
In the comparison Control vs CV + R, the MDS plot also demonstrates a good discrimination between the two experimental groups, with Control located in the region corresponding to higher values of leading logFC dimension 1 and CV + R with lower values of leading logFC dimension 1 and dimension 2 (Fig. 3A). This information is corroborated by the correlation coefficients presented in Fig. 3C. Some of the identified DEGs are presented as red dots in the volcano plot (Fig. 3B).
(A) Multi-dimensional scaling (MDS) plot showing the relation between Control and CV + R diet samples. (B) Volcano plot showing the relationship between FC and evidence of differential expression (-log FDR) for Control vs CV + R. (C) Clustering of sample to sample correlations—Control vs CV + R. Hierarchical clustering of treatment condition samples based on Pearson (left) and Spearman correlations (right).
In the comparison Control vs CV + M, less discrimination was observed between these two experimental groups in the MDS plot (Fig. 4A) and in the correlations (Fig. 4C), when comparing with previous comparisons. Thus, a much lower number of DEGs was identified between these two experimental groups, marked as red dots in the volcano plot (Fig. 4B) compared with Control vs CV and Control vs CV + R.
(A) Multi-dimensional scaling (MDS) plot showing the relation between Control and CV + M diet samples. (B) Volcano plot showing the relationship between FC and evidence of differential expression (-log FDR) for Control vs CV + M. (C) Clustering of sample to sample correlations—Control vs CV + M. Hierarchical clustering of treatment condition samples based on Pearson (left) and Spearman correlations (right).
Functional classification of DEGs
In order to understand which biological pathways and processes were represented by the identified DEGs, a functional analysis was performed. For this analysis, in comparisons Control vs CV and Control vs CV + R, were selected the 50 most up and 50 most downregulated genes and were also included the genes involved in fatty acid metabolism and oxidative stress. This analysis was not possible to perform in Control vs CV + M due to the reduced number of identified DEGs. Table 1 displays the identified significant DEGs for comparisons Control vs CV and Control vs CV + R, its respective functional description, and Log2 FC.
Table 1 Identified DEGs by functional analysis with Cytoscape software for comparisons Control vs CV and Control vs CV + R, its respective functional description and Log2 FC.
In Control vs CV (Fig. 5A), the majority of the DEGs are downregulated in the Control group and consequently, upregulated in CV and are marked as blue circles. Only the following DEGs: SCD, BRCA1, CES1, MGST2, APOD, LPL, C3, and HADHA are upregulated in Control and are marked as red circles. In this comparison, DEGs are associated to the following biological processes (squares): “Fatty acid catabolic process” (P = 1.3E-08), “Fatty acid biosynthetic process” (P = 4.9E−13), “Cholesterol metabolic process” (P = 5.4E-04), “Negative regulation of lipid metabolic process (P = 1.9E-06), “Triglyceride metabolic process (P = 8.8E-05), and “Antioxidant activity” (P = 0.008), and are associated to the following pathways (diamonds): “Mitochondrial fatty acid β-oxidation” (P = 2.4E-08), and “Glycerophospholipid biosynthesis” (P = 0.002).
Functional analysis performed with Cytoscape software using the 50 most up and 50 most downregulated genes and genes of interest for (A): Control vs CV comparison and (B): Control vs CV + R comparison. The significant DEGs and the interactions with their related pathways and biological processes are presented in the figure. Legend: squares = biological processes; diamonds = reactome pathways; circles = DEGs; shape size = according to the p-value of the term in its own group; red colour = upregulated in Control; blue colour = downregulated in Control.
Regarding Control vs CV + R (Fig. 5B), it was observed that all DEGs associated with “Peptide chain elongation” (P = 0.03) and “Selenocysteine synthesis “ (P = 0.02) pathways are upregulated in Control (red diamonds). For all other identified biological processes and pathways, the majority of the DEGs are downregulated in the Control group (blue circles), with the exception of SCD and NLRP3 (red circles). The remaining biological processes (squares) identified are: “Protein transmembrane import into intracellular organelle” (P = 0.02), “Regulation of fatty acid metabolic process” (P = 1.8E-05), “Regulation of cholesterol metabolic process” (P = 6.9E-04), “Superoxide dismutase activity” (P = 4.6E-04), and “Cellular response to superoxide” (P = 0.009). The identified pathways (diamonds) are: “Fatty acid metabolism” (P = 3.3E−12), “Peroxisomal lipid metabolism” (P = 6.6E-06), “Protein localization” (P = 3.6E-04), “Peroxisomal protein import” (P = 5.4E-07), and “Mitochondrial fatty acid β-oxidation” (P = 2.3E-05).
In the comparison Control vs CV + M, three DEGs were identified. The ncRNA gene MIR10390 and the protein-coding gene LMO7 are upregulated in Control, while the protein-coding gene IGFN1 is upregulated in the CV + M group.
A total of 1329 protein identifications were obtained in the muscle of finishing pigs following the LC–MS analysis (Supplementary file 2). The following sections describe the differentially abundant proteins between experimental groups (CV, CV + R, CV + M) and control.
Control vs CV
21 proteins were differentially abundant between CV and control groups (Supplementary Table S2). The former had 11 proteins with higher abundance. These proteins included acyl-CoA binding domain containing 7 (ACBD7), and ribosomal protein S25 (RPS25), with fatty-acyl-CoA binding (GO:0000062), and structural constituent of ribosome (GO:0003735) molecular functions. It also had a higher abundance of contractile apparatus proteins including troponin T (tnnt3), and myosin light-chain 1 (MYL1). Conversely, control had higher abundance of other structural muscle proteins, myosin-2 (MYH2), and ACTN1. This group had also higher abundance of mRNA metabolism proteins, including 2-iminopropanoate deaminase (RIDA), and ribonucleoprotein D (HNRNPD) with mRNA catabolic process (GO:0006402), and positive regulation of translation (GO:0045727) biological process annotations, respectively.
Control vs CV + R
This comparison has yielded 18 differentially abundant proteins, with 5 of them being highly abundant in CV + R (Supplementary Table S3). This latter group of proteins had heterogeneous biological process annotations including respiratory electron transport chain (GO:0022904, Rieske domain-containing protein -UQCRFS1), and fatty acid β-oxidation (GO:0006635, hydroxysteroid 17-β dehydrogenase 10—HSD17B10). In turn, control increased the abundance of several myosin proteins (MYH1, MYH2, MYH4 and MYH7), which are involved in muscle contraction/motor activity (GO:0003774). It also increased the abundance of ATP binding (GO:0006883) proteins such as sodium/potassium-transporting ATPase subunit alpha (ATP1A3), and creatine kinase (CKMT1), with sodium/potassium homeostasis (GO:0006883, GO:0030007), and phosphocreatine biosynthetic process (GO:0046314) biological process annotations, respectively.
Control vs CV + M
This comparison had the lowest number of differentially abundant proteins, with 5 and 7 highly abundant proteins in CV + M and control, respectively (Supplementary Table S4). The former group had higher abundance of ribosome (RPS17), and SERPIN domain-containing (SERPINA3-2) proteins which participate in translation (GO:0006412), and negative regulation of endopeptidase activity (GO:0010951), respectively. It also increased the abundance of SPARC, a protein that regulates cell growth (anatomical structure development, GO:0048856). Control increased the abundance of myosin-binding protein H (MYBPH), and ACTN1, both being contractile apparatus proteins. In addition, they also increased the abundance of NADH dehydrogenase ubiquinone iron-sulphur protein 3 (NDUFS3), which partakes in mitochondrial electron transport (GO:0006120), and creatine kinase (CKMT1).
A supervised analysis was conducted using the DIABLO procedure to detect multi-omics signatures that distinguish the experimental groups. A threshold of 0.7 was used to illustrate significant correlations between genes and proteins, as seen in the circosplots of Figure S3.
Control vs CV
In this comparison, the two datasets presented a high correlation of 0.94 (Figure S4A). The correlation circle plot obtained for this comparison demonstrates three separate clusters of genes and proteins that positively correlate with each other (Figure S5A). These are depicted in the circosplot obtained in Figure S3A, which shows a high number of positive correlations. Two contractile apparatus proteins, ACTN1 (I3LLY3) and MYBPH (I3LIE7), actin and myosin-binding proteins, are positively correlated with genes including NCAPD3 (participates in cell division), SCD (participates in lipid metabolism), and ZNF879 (regulates transcription).
Histone H4 (P62802) is positively correlated with FANCD2, and NCAPD3, two genes with maintenance of chromosomal stability and histone binding functions, respectively.
The ribosomal protein RPS25 (F2Z5G8), and the collagen-binding SPARC protein (P20112) were negatively correlated with the PAPPA2 gene. The latter is involved in proteolysis and regulates cell growth.
Control vs CV + R
This comparison yielded a high dataset correlation of 0.95, similarly to the CV vs control comparison (Figure S4B). As seen in Figure S5B and Figure S3B, there is a negative correlation of a cluster of genes with a particular protein, CAMK2A (I3LNG5). These genes include BAX (involved in apoptosis), DDA1 (positively regulates protein catabolism), and MAD2L2 (negatively regulates transcription). The aforementioned protein is a kinase that activates transcription in developing neurons. Interestingly, another protein (HNRNPD—F1RVC9), which regulates gene expression, is positively correlated with these genes. Both of these proteins are highly abundant in control.
In turn, three myosin proteins (MYH1, MYH4, MYH7), NEB (actin binding), and CRYAB (muscle development) are positively correlated with lipid metabolism genes such as CIDEC (lipid droplet organisation/biogenesis), and SCD (phospholipid, cholesterol and triglyceride synthesis). These genes and proteins were all highly abundant in the control group compared to CV + R.
CV + M vs control
This comparison yielded the lowest number of differentially expressed genes and abundant proteins, achieving the lowest correlation between the two datasets (Figure S4C). As seen in the correlation circle plot (Figure S5C) and circosplot (Figure S3C), there is a positive correlation between the IGFN1 gene and three proteins: SPARC (P20112), SERPIN domain-containing protein (F1SCC6), and TMEM263 (A0A287B182) proteins. Likewise, there is another positive correlation between two genes (LMO7 and mir10390) and four proteins: ACTN1 (I3LLY3), NDUFS3 (A0A286ZNN4), CAMK2A (I3LNG5), and TUBA4A (F2Z5S8).