Identification of hub genes in meat quality of grass carp (Ctenopharyngodon idellus) fed with faba bean by muscle tissue transcriptomic analysis

Abstract

Crisp grass carp is a type of grass carp (Ctenopharyngodon idella) that has been fed faba beans. Flesh quality and expression of muscle quality-related genes in muscle of grass carp were evaluated. Compared with ordinary grass carp, 1901 differentially expressed genes (DEGs) were up-regulated and 1628 DEGs were down-regulated in crisp grass carp. The GO analysis results showed that the DEGs were enriched in ribosomes and proteins. Kyoto Encyclopaedia of genes and genomes (KEGG) pathway analysis results showed that amino acid metabolism was closely related to meat quality. COX3, ACADS, ACAT1, ECHS1, EHHADH, HADHAA, HADH, HMGCL, MCCC1 and MCCC2 were the top 10 hub genes most closely associated with other nodes and showed a strong positive correlation. Hardness was positively correlated with the COX3, ECHS1, ACADS, ACAT1, HADHA, HADH, HMGCL and MCCC1. These results have been demonstrated to play important roles in clarifying the molecular mechanism of muscle firmness increase in crisp grass carp.

Highlights

• Feeding faba bean meal improved texture parameters of grass carp.

• A more comprehensive crisp grass carp transcriptome database was established.

• Several amino acid metabolism pathway that were closely related to meat quality were identified.

Keywords:

• Ctenopharyngodon idellus
• faba bean
• meat quality
• RNA sequencing

Introduction

Grass carp (Ctenopharyngodon idellus) is a kind of freshwater fish which is widely farmed in China and has great economic importance (Yang et al. 2020). Grass carp is well known for excellent adaptability, rapid growth and large size, and it is one of the ‘Four Domesticated Fish’ in China (Liu et al. 2013). Through digestion and absorption, it can convert plant protein into animal protein directly (Lu et al. 2020). Grass carp is one of the most productive freshwater fish in China, which is also the species fish with the highest yield in freshwater fish farming worldwide (Rao and Su 2015; Wu et al. 2022).

Fish is an essential portion of the human diet because it is an excellent source of protein and unsaturated fatty acids, both of which are thought to be beneficial to human health. While breeding output is rising, it frequently appears that meat quality is declining, which has a significant negative impact on the meat’s economic value. Numerous studies have been attempted to enhance the quality of the meat from farmed grass carp (Yu et al. 2020; Yang et al. 2022; Li et al. 2023). Among them, feeding grass carp with faba bean to embrittle its meat quality is an innovative method that enhances the meat quality of grass carp (Zhou et al. 2019). Grass carp raised by feeding faba bean are called crisp grass carp (Lin et al. 2009). Studies have shown that the muscle characteristics of crisp grass carp have greatly changed compared with ordinary grass carp, in which the hardness, chewiness and springiness have increased significantly, the diameter of muscle fibres has decreased, and the density of muscle fibres has increased (Lin et al. 2009), which has increased its added value. Faba beans are also known to increase muscle hardness in fish such as Dicentrarchus labrax (Adamidou et al. 2009) and Oreochromis niloticus (Li et al. 2023). Fava beans also affect the fatty acid content of fish, increasing fat deposits in the viscera (Tian et al. 2019).

In the previous studies of crisp grass carp muscle, it has been found that the chewiness, flesh hardness, adhesiveness and elasticity of grass carp were significantly increased by feeding the fish faba beans (Yang et al. 2015), and the amino acid and fatty acid composition were influenced by fava bean in grass carp (Jiang et al. 2020; Li et al. 2020; Xu et al. 2020). The amount of collagen was another factor affecting the firmness of fish meat. A higher content of collagen produces a firmer raw fish meat, and the increase in the expression of type I collagen in crisp grass carp is higher than those of grass carp (Yu et al. 2014, 2019). However, long-term intake of whole fava bean can lead to poor palatability of crisp grass carp, oxidative stress, intestinal inflammation and decreased growth rate (Yu et al. 2014, 2017; Ma et al. 2020), and the maximal proportion of dietary protein replacement with faba bean meal was 70% based on the flesh quality of grass carp (Fu et al. 2022). Therefore, it is expedient to analyse systematically muscle firmness increase of crisp grass carp in the gene levels. Understanding the impact of faba beans on skeletal muscle quality is important for optimising feed used for high-quality crisp grass carp production.

Improving the muscle mass of fish is a major goal of aquaculture. Transcriptome analysis can accurately, efficiently and quickly obtain transcriptome sequence information and expression of specific cells or tissues in a certain state, as well as analyse changes in gene expression (Azaman et al. 2019; Gill and Dhillon 2022), which provides an effective way to explore the molecular mechanism of faba bean regulating grass carp meat quality. This study intends to perform transcriptome sequencing on the muscle tissue of crisp grass carp and ordinary grass carp and to screen out differentially expressed genes (DEGs) and signalling pathways related to muscle, so as to lay a foundation for analysing the molecular mechanism of faba bean regulating the development of grass carp meat quality.

Materials and methods

Tissue collection of experimental fish

Healthy grass carp were purchased from an aquaculture farm in Zunyi, Guizhou Province, China. The fish were first temporarily cultured in a cement pond (5 m × 5 m × 1.5 m) for 1 week and the feed amount for each day was 2–3% of fish weight. A total of 180 fish with initial weight of 768 ± 75 g were randomly divided into crisp grass carp (group A) and ordinary grass carp groups (group B), with three replicates each group. They were cultured in six cement ponds (2 × 2 × 1.5 m), with 30 fish in each pond. Crisp grass carp were fed solely with whole faba beans, the feed was soaked in about 0.15% salt water for 24 h, and then soaked in water for 12 h until the broad beans were opened after the germ. The single feeding amount of broad bean accounted for 2–3% of the body weight of fish. The ordinary grass carp were fed with commercial diet (crude protein: 329.9 g kg⁻¹; crude lipid: 43.8 g kg⁻¹; Tongwei, Chengdu, China). The fish were fed twice per day (at 8:00 and 17:00). The water temperature was kept at 25–30 °C, pH was 6.5–7.5, and dissolved oxygen was above 5.0 mg/L. The final weights of crisp grass carp and ordinary grass carp were 1992 ± 125 g and 2457 ± 132 g after 120 days, respectively. One fish was randomly selected from each pond to collect its muscle tissue, which was placed at −80 °C until RNA extraction. The procedure in this experiment was approved by the Ethics Committee of Experimental Animal of Zunyi Normal College (protocol number: ZUNSHIFA[2018]08).

Determination of muscle physical parameters and texture characteristics

Five fish from each group were used for determination of muscle physical and texture parameters. The muscle pH was analysed following the method reported by Knowles et al. (2007). Cooking loss was measured as follows: 2–3 g of dorsal muscle samples were dried with filter papers and weighed (weight was W₀), and then heated for 1 min in a pan of boiling water. The samples were weighed again (weight was W₁) after cooling to room temperature and removing surface water. The calculation formula was as follows: cooking loss (%) = 100 × (W₀ – W₁)/W₀. The texture of the dorsal muscle, including hardness, springiness, chewiness and cohesiveness, was performed with the texture profiles analysis (TPA) using a Universal TA device (Tengba, Shanghai, China). Test conditions were as follows: a 25 mm × 25 mm flat-bottomed cylindrical probe, a compression ratio of 30%, a test speed of 1 mm/s and a post-test speed of 1 mm/s with the staying time of 2 s. Test conditions of shear force were as follows: compression distance of 20 mm, a test speed of 1 mm/s and a post-test speed of 1 mm/s with a staying time of 2 s.

RNA extraction and sequencing

Total RNA was isolated using the Trizol Reagent (Invitrogen Life Technologies, Carlsbad, CA), after which the concentration, quality and integrity were determined using a NanoDrop spectrophotometer (Thermo Scientific, Waltham, MA). Three micrograms of RNA were used as input material for the RNA sample preparations. The mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. The cDNA was synthesised using the mRNA fragment as templates. To select cDNA fragments of the preferred 400–500 bp in length, the library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, CA). DNA fragments with ligated adaptor molecules on both ends were selectively enriched using Illumina PCR Primer Cocktail in a 15 cycle PCR reaction. Products were purified (AMPure XP system, Beckman Coulter, Beverly, CA) and quantified using the Agilent high sensitivity DNA assay on a Bioanalyser 2100 system (Agilent, Santa Clara, CA). The sequencing library was then sequenced on NovaSeq 6000 platform (Illumina, San Diego, CA).

Transcriptome assembly and reads mapping

Samples are sequenced on the platform to get image files, which are transformed by the software of the sequencing platform, and the original data in FASTQ format (Raw Data) are generated. Sequencing data contain a number of connectors, low-quality Reads, so we use Cutadapt (v1.15) software to filter the sequencing data to get high quality sequence (Clean Data) for further analysis. The reference genome and gene annotation files were downloaded from genome website. The filtered reads were mapping to the reference genome using HISAT2 v2.0.5.

Differential expression analysis

HTSeq (0.9.1) statistics was used to compare the Read Count values on each gene as the original expression of the gene, and then used FPKM to standardise the expression. Then difference expression of genes was analysed by DESeq (1.30.0) with screened conditions as follows: expression |log2FoldChange| >1, significant p value <.05. At the same time, we used R language Pheatmap (1.0.8) software package to perform bi-directional clustering analysis of all different genes of samples. We got the heatmap according to the expression level of the same gene in different samples and the expression patterns of different genes in the same sample with Euclidean method to calculate the distance and Complete Linkage method to cluster.

GO and KEGG enrichment analysis

All the genes to terms in the Gene Ontology database were mapped and the numbers of differentially enriched genes in each term were calculated. Using topGO to perform GO enrichment analysis on the differential genes, calculate p value by hypergeometric distribution method (the standard of significant enrichment is p value <.05), and find the GO term with significantly enriched differential genes to determine the main biological functions performed by differential genes. ClusterProfiler (3.4.4) software was used to carry out the enrichment analysis of the Kyoto Encyclopaedia of genes and genomes (KEGG) pathway of differential genes, focusing on the significant enrichment pathway with p value <.05.

Interaction analysis of differential gene protein network

The STRING database (https://string-db.org/) is used for protein interaction analysis to reveal the relationship between target genes. The protein–protein interaction (PPI) network of the DEGs was constructed based on the Search Tool for the Retrieval of Interacting Genes (STRING) database22 (STRING, version 9.1, http://string91.embl.de/), and predicted by the use of Cytoscape23 (version 3.0; http://cytoscape.org/). The hub genes were then selected with a connectivity degree 10 after calculating the degree of each node. Module analysis of the PPI network was performed with the parameters of minimum size >4 and minimum density <0.05 using ClusterONE.

Correlation analysis and qRT-PCR validation

To validate the results of the transcriptome data, top 10 hub genes were randomly selected and their mRNA levels were measured by qRT-PCR. Specific primer sequences for qRT-PCR were as follows: peroxisomal acyl-coenzyme A oxidase 3 (COX3, 5F: CCAGTCCCTCACCCTAACAA; R: GTGGTGTTCGGAGGTGAAGT), enoyl-CoA hydratase (ECHS1, 5F: TAAGGCCTTGGACACCTTTG; R: TTTGTTCCAGTGTGCCAAAA), short-chain specific acyl-CoA dehydrogenase (ACADS, 5F: TCGACAGGAGTCATCGTCAG; R: TTGAGAACCCATTCGTTTCC), acetyl-CoA acetyltransferase (ACAT1, 5F: TGTCATCGTCAGTGCTGTCA; R: AGCTTGTCTGGTTGGTGCTT), peroxisomal bifunctional enzyme (EHHADH, 5F: GGGTGTGTAGACCAGCCACT; R: ATCCTCTTCCCAAGGCTCAT), trifunctional enzyme subunit alpha (HADHA, 5F: TGCCATCAAGAGTGCTGTTC; R: GGCTACAATGGGTTTGGATG), trifunctional enzyme subunit alpha (HADH, 5F: TAAATCTGCGAAGGGGATTG; R: ACAATAGCCTCCACCACCAG), hydroxymethylglutaryl-CoA lyase (HMGCL, 5F: CTCGGTGGAGGAGAGTCTTG; R: TGTCACCCAATGAAACCTCA), methylcrotonyl-CoA carboxylase subunit alpha (MCCC1, 5F: CAACATCGACTTCCTGCTGA; R: TCGCTGGTCTGCTCTCTGTA) and methylcrotonyl-CoA carboxylase beta chain (MCCC2, 5F: AGCCAGAGTTCGTCCTTTGA; R: GTGCTTGCATTCGCTCATAA). The total mRNA of the samples was extracted using the Trizol Kit. PrimeScript™ RT Reagent Kit with gDNA Eraser (Takara, Kusatsu, Japan) was then used for cDNA synthesis based on manuals. qRT-PCR was conducted with a Roche Step One real-time PCR system (Roche, Indianapolis, IN). PCR amplification of all samples was performed in triplicate. β-Actin (5F: TGCCATGTATGTGGCCATCC; R: TCTTTCGGCTGTGGTGGTGA) was used as the normalisation control to calculate the relative expression levels (Li et al. 2019). The melting curve was analysed to confirm the target specificity. The 2^−ΔΔCT method was used to calculate the relative gene expression values.

Results

Muscle physical parameters and texture characteristics

In this study, the muscle physical parameters and texture characteristics for ordinary grass carp and crisp grass carp are given in Table 1. Lower water holding capacity was observed in crisp grass carp group when compared with ordinary grass carp group (p < .05). Hardness, springiness and gumminess in crisp grass carp group were significantly (p < .05) higher compared with ordinary grass carp. Besides, muscle pH, cooking loss, chewiness and cohesiveness showed no obvious difference between the two groups (p > .05).

Table 1. Muscle physical and texture parameters of grass carp after 120-days feeding experiment.

Parameters	Crisp grass carp	Ordinary grass carp
pH	6.56 ± 0.09	6.66 ± 0.01
Water holding capacity, %	7.19 ± 1.62^b	10.94 ± 2.49^a
Cooking loss, %	24.12 ± 2.33	27.03 ± 2.58
Hardness, kg	1.23 ± 0.29^a	0.75 ± 0.12^b
Springiness, mm	4.98 ± 0.01^a	3.08 ± 0.56^b
Gumminess, kg	0.89 ± 0.08^a	0.61 ± 0.15^b
Chewiness, mJ	38.18 ± 5.98	29.56 ± 6.19
Cohesiveness, mm	0.55 ± 0.03	0.49 ± 0.07

Results are presented as the mean ± SEM (n = 5). Values with different letters in the same line indicate significant differences (p < .05).

Transcriptome sequencing and assembly results

In this study, a total of 40,564,018–44,425,106 raw reads and 37,910,970–41,601,388 clean reads were obtained through transcriptome analysis of six muscle tissue samples. The alignment rate of clean reads from ordinary grass crap and crisp grass crap mapped to the reference genome ranged from 94.28% to 95.14%. The Q20 is higher than 98.06% and the Q30 is higher than 94.49% (Table 2). The generated sequence data have been submitted to the NCBI SRA database (accession number: PRJNA1053680).

Table 2. Summary statistics of transcriptome sequencing from muscle tissue between ordinary grass carp and the crisp carp group.

Sample	Raw reads	Raw bases, bp	Clean reads	Clean data, bp	Total_Mapped	Q20, %	Q30, %
Group_A1	40,564,018	6,125,166,718	37,910,970	5,724,556,470	94.82%	98.06	94.49
Group_A2	44,425,106	6,708,191,006	41,601,388	6,281,809,588	94.82%	98.41	95.26
Group_A3	42,834,600	6,468,024,600	39,930,392	6,029,489,192	94.77%	98.22	94.83
Group_B1	44,106,208	6,660,037,408	41,239,048	6,227,09,6248	94.91%	98.35	95.19
Group_B2	44,332,974	6,694,279,074	41,533,704	6,271,589,304	95.14%	98.16	94.7
Group_B3	43,989,620	6,642,432,620	41,083,156	6,203,556,556	94.28%	98.26	94.93

Group_A1, Group_A2 and Group_A3: Crisp carp group; Group_B1,Group_B2 and Group_B3: Ordinary grass carp.

Differentially expressed gene identification and analysis

Gene expression was calculated according to the FPKM algorithm. Differential analysis of gene expression was performed using DESeq, and screening was performed according to |log2FoldChange| >1 and p value < .05. There were 3529 DEGs in crisp grass carp (group A) compared with ordinary grass carp groups (group B). Among them, 1901 differential genes were down-regulated, and 1628 differential genes were up-regulated (Figure 1). In the DEGs of the top 20 up-regulated, large neutral amino acids transporter small subunit 4 (SLC43A2) is the gene with the most significant difference. In the DEGs of the top 20 down-regulated, l-lactate dehydrogenase A (LDHA) is the gene with the most significant difference. Among all DEGs, some genes related to muscle growth and development were also screened out, such as insulin-like growth factor-binding protein 2, myosin heavy chain (Table 3). Differentially expressed genes (DEGs) identified in the muscle tissue between ordinary grass carp and the crisp carp group (Top20) are presented in Figure 2. GO enrichment analysis showed that DEGs were mainly enriched in collagen type I trimer, l-amino acid transmembrane transporter activity, carboxylic acid transmembrane transporter activity, amino acid transmembrane transporter activity, l-leucine transmembrane transporter activity, neutral amino acid transmembrane transporter activity, neutral amino acid transport, amino acid transport, l-alpha-amino acid transmembrane transport, branched-chain amino acid transmembrane transporter activity, leucine import across plasma membrane, amino acid transmembrane transport and collagen metabolic process. KEGG enrichment analysis showed that DEGs were mainly enriched in protein digestion and absorption, cysteine and methionine metabolism, glycolysis/gluconeogenesis, biosynthesis of amino acids, adipocytokine signalling pathway, tryptophan metabolism and pyruvate metabolism.

Figure 1. Volcano map and heat map of clustering of differentially expressed genes in Ctenopharyngodon idellus.

Figure 2. Classify and analysis according to gene function. (A) GO enrichment of top 20 up-regulated and down-regulated DEGs. (B) KEGG enrichment of top 20 up-regulated and down-regulated DEGs.

Table 3. Differentially expressed genes (DEGs) identified in the muscle tissue between ordinary grass carp and the crisp carp group (Top20).

Gene name	log2FoldChange	pval	padj	Up/down
SLC43A2	5.064656592	6.80492E–33	1.08032E–28	Up
SLC15A4	4.705924222	9.08516E–33	1.08032E–28	Up
CTSL	4.48421604	2.04893E–32	1.62425E–28	Up
ASB18	4.302901408	9.31283E–27	5.53694E–23	Up
PRODH	4.288735951	3.07603E–24	1.21923E–20	Up
RETREG1	3.646621252	8.11353E–24	2.75651E–20	Up
SYT13	3.444645426	6.42403E–22	1.69751E–18	Up
TMEM117	5.606150055	8.53125E–21	1.44922E–17	Up
KLHL38	3.684010648	1.02518E–20	1.62538E–17	Up
ALSIN2	3.786370075	2.1948E–20	3.2623E–17	Up
FBXO30	4.343309036	3.60615E–20	4.76452E–17	Up
SLC16A9	3.556296732	4.43136E–20	5.54666E–17	Up
CCDC129	4.02000405	2.2623E–19	2.562E–16	Up
MSTO1	3.413871711	7.10594E–19	7.34755E–16	Up
FBXO32	3.077803356	7.81533E–19	7.74434E–16	Up
SLC7A8	6.709140063	3.72444E–18	3.54299E–15	Up
NFKBIE	3.46188969	1.66097E–17	1.46301E–14	Up
ABHD2B	3.200719079	4.2128E–17	3.45478E–14	Up
SLC1A5	2.850515434	7.89862E–17	6.05951E–14	Up
ARRDC2	2.88068052	1.1498E–16	8.54514E–14	Up
LDHA	−4.36656894	1.39724E–26	6.64581E–23	Down
AGL	−3.66368887	8.47661E–23	2.51988E–19	Down
SMARCD3	−3.80962945	7.67351E–22	1.82491E–18	Down
ABI3BP	−3.708319297	9.88402E–22	2.13693E–18	Down
TP53BP1	−3.619697977	7.4131E–21	1.38612E–17	Down
MR1	−8.336943833	7.57699E–21	1.38612E–17	Down
COL1A1	−3.294263698	3.21588E–20	4.49883E–17	Down
PDK3	−3.414423635	6.83228E–20	8.12427E–17	Down
FNDC1	−3.225926906	2.96302E–17	2.51666E–14	Down
AAMM50B	−4.606184319	7.56746E–17	5.99898E–14	Down
COL1A2	−2.887803272	2.56396E–16	1.7437E–13	Down
PGK1	−2.819612167	5.67971E–16	3.57571E–13	Down
TNNT3	−2.865108634	7.9281E–16	4.71365E–13	Down
HP1BP3	−3.648142909	1.11962E–15	6.33972E–13	Down
AHSA1	−3.432529968	7.74482E–15	4.09305E–12	Down
ADH1	−3.698975919	9.94866E–15	5.14346E–12	Down
BCL2L13	−2.932306208	1.93091E–14	9.62513E–12	Down
SMTN	−2.671925777	3.06599E–14	1.45831E–11	Down
KYAT1	−3.349710817	4.46425E–14	2.08174E–11	Down
MYOM1	−2.728438474	6.16919E–14	2.71696E–11	Down

GO enrichment analysis of DEGs

The results of GO enrichment analysis showed that DEGs were mainly enriched in 5484 GO terms of biological process (BP), molecular function (MF) and cellular component (CC), among which BP was the most prevalent (3493). In the BP, the DEGs are mainly involved in BP and cellular process, containing 1072 and 916 DEGs, respectively. In the MF, the DEGs involved in the MF process are the most, including 1030 DEGs. In the CC, the DEGs are mainly involved in the process of CC and cellular anatomical entity, including 967 and 927 DEGs, respectively.

GO enrichment analysis showed that the top 20 terms enriched DEGs were enriched in proteins and ribosomes (Figure 3(A)). GO enrichment showed that up-regulated genes were significantly enriched in the ribosomes, including ribonucleoprotein complex (GO:1990904), large ribosomal subunit (GO:0015934), small ribosomal subunit (GO:0015935), cytosolic small ribosomal subunit (GO:0022627), cytosolic large ribosomal subunit (GO:0022625), ribosomal subunit (GO:0044391), cytosolic ribosome (GO:0022626), ribosome (GO:0005840) and structural constituent of ribosome (GO:0003735). Protein-related GO terms include peptide biosynthetic process (GO:0043043), peptide metabolic process (GO:0006518) and translation (GO:0006412) (Figure 3(B)).

Figure 3. (A) GO enrichment of DEGs. (B) GO enrichment of up-regulated DEGs. (C) GO enrichment of down-regulated DEGs.

GO enrichment analysis showed that down-regulated genes were significantly enriched in the protein catabolism, including proteasome accessory complex (GO:0022624), proteolysis (GO:0006508), modification-dependent protein catabolic process (GO:0019941), ubiquitin-dependent protein catabolic process (GO:0006511), proteasome regulatory particle (GO:0005838), organonitrogen compound catabolic (GO:1901565), process cellular protein catabolic process (GO:0044257), proteolysis involved in cellular protein catabolic (GO:0051603), protein catabolic process (GO:0030163), peptidase complex (GO:1905368), endopeptidase complex (GO:1905369) and proteasome complex (GO:0000502) (Figure 3(C)).

KEGG enrichment analysis of DEGs

KEGG enrichment analysis showed that DEGs in muscle tissue were enriched into 252 known KEGG pathways, of which 45 signalling pathways were significantly enriched (p < .05).

The top 20 KEGG pathways enriched by DEGs are shown in Figure 4(A), and nine of the 20 top pathways were related to metabolism of energy, protein, amino acid and fatty acid. The top 20 KEGG pathways enriched by up-regulated differential genes are shown in Figure 4(B). There were six KEGG pathways related to amino acid metabolism, including arginine and proline metabolism, beta-alanine metabolism, lysine degradation, glutathione metabolism, tryptophan metabolism and valine, leucine and isoleucine degradation. Ribosome is also an upregulated GO term. The top 20 KEGG pathways enriched by down-regulated differential genes are shown in Figure 4(C). There were two pathways related to protein, including protein processing in endoplasmic reticulum and ubiquitin mediated proteolysis.

Figure 4. (A) KEGG pathway enrichment of DEGs. (B) KEGG pathway enrichment of up-regulated DEGs. (C) KEGG pathway enrichment of down-regulated DEGs.

PPI network construction

The PPI network with 72 nodes was constructed by analysing DEGs associated to differential amino acid metabolic pathways using the STRING database (Figure 5(A)). COX3, ECHS1, ACADS, ACAT1, EHHADH, HADHA, HADH, HMGCL, MCCC1 and MCCC2 were the top 10 hub genes most closely associated with other nodes, according to further research with cytoHubba (Figure 5(B)).

Figure 5. PPI network construction and hub gene screening.

Correlation analysis and qRT-PCR validation

The correlations of top 10 hub genes with meat quality are shown in Figure 6(A). Hardness was positively correlated with the COX3, ECHS1, ACADS, ACAT1, HADHA, HADH, HMGCL and MCCC1. There is a strong positive correlation among the top 10 hub genes. To validate the RNA-Seq results obtained, 10 genes were randomly selected for validation by qRT-PCR. These genes that were involved in the amino acid metabolism were closely related to meat quality. As shown in Figure 6(B), the qPCR expression patterns of the selected genes agreed with the results of RNA-Seq analysis, indicating the reliability and accuracy of the RNA-Seq method used in the present study.

Figure 6. (A) Correlation analysis between the top 10 hub genes and meat quality. (B) Comparison between RNA-Seq result and qRT-PCR validation result. The y-axis is the relative fold change. Data are presented as mean ± SEM.

Discussion

Changes in muscle texture are generally regarded as an important indicator of the ‘embrittlement’ of crisp grass carp, especially the chewiness of muscles, which are an important factor for the ‘crispy’ of muscles. Feeding faba beans is an effective way to improve the meat quality of grass carp. Our results also showed that some texture parameters of crisp grass carp, such as the hardness, springiness and gumminess were significantly improved, while the water holding capacity in muscle was reduced. These results indicated that feeding faba bean was an effective way for improving muscle texture of grass carp.

SLC43A2 mediates the transport of leucine, phenylalanine, valine and methionine, which is critical for early development in animal (Guetg et al. 2015). LDHA can catalyse the conversion of l-lactic acid into pyruvate, which is closely related to the generation of ATP (Feng et al. 2018). Myosin, which contains components like MHC and MLC, is one of the primary structural and functional proteins that make up skeletal muscle and cardiac muscle (Agarwal et al. 2020). The size, quantity and type of muscle fibres are all influenced by MHC during the development of skeletal muscle. In the present study, the expression of 1628 genes in crispy grass carp was significantly up-regulated, while the expression of 1901 genes was significantly down-regulated. In the DEGs of the top 20 up-regulated, SLC43A2 is the most significant difference. In the DEGs of the top 20 down-regulated, LDHA is the most significant difference. Genes of type I collagen and type II collagen were once considered tightly related to muscle-hardening in grass carp fed with faba bean (Yu et al. 2019), and we have identified these genes in the current study. These DEGs may be related to the muscle growth and development of grass carp, which is worthy of further digging and research.

In this study, eight of the top 20 up-regulated GO terms were related to ribosomes. In fact, ribosomes are the main place for cells to synthesise proteins, and the process of RNA translation into proteins is completed in ribosomes (Zimmermann et al. 2016), whereas eukaryotic 80S ribosomes are composed of 40S small and 60S large subunits. The small subunit decodes transcripts by promoting base-pairing between codons within an mRNA and their corresponding transfer RNAs (tRNAs) (Buszczak 2023). The large subunit contains the peptidyl transferase centre which catalyses peptide bond formation between amino acids (Buszczak 2023). Most of the protein in skeletal muscle consists of the contractile proteins actin and myosin. The myocyte hypertrophy is independent of the proliferation of satellite cells and myonuclear accretion, in contrast to long-term remodelling and maintenance of skeletal muscle homeostasis, which appear to require satellite cells together with their influence on non-skeletal muscle cells (Fry et al. 2017). Our results showed that ribosome was the most DEGs among the top 20 up-regulated KEGG pathways, suggesting that the muscle cell coordinates its anabolic and nonanabolic processes in a complex and sophisticated manner, and ribosome biogenesis is important for hypertrophic growth of skeletal muscle. At the same time, translation and transcription are clearly linked and preferably synchronised during this process, probably to optimise the response to energy demands and the potential for growth and muscle adaptation (von Walden 2019). In addition, ribosomal proteins have been shown to influence cell differentiation, proliferation and death as well as play an essential component in muscle growth (Genuth and Barna 2018), and our results suggest that ribosomal proteins play an important role in maintaining the body’s normal life activities, and the content of amino acids is related to muscle nutrition and quality.

Meat quality, including amino acid content, is an important factor affecting meat flavour and consumer acceptance (Zhang et al. 2022). Amino acid concentration is an essential indicator for assessing the quality of meat and is crucial in enhancing meat flavour (Huo et al. 2021; San et al. 2021). It has been reported that in farmed fish, the composition of tissue amino acid is greatly affected by dietary nutrients. This is consistent with those previous studies, among the up-regulated pathways caused by feeding faba beans, there are six KEGG pathways related to amino acid metabolism, indicating that arginine, proline, β-alanine, lysine, tryptophan, leucine and isoleucine acid play an important role in improving grass carp meat quality. Some KEGG analysis results of this study are similar to those of other fishes, such as glycolysis/gluconeogenesis (Li et al. 2017), arginine and proline metabolism (Canosa and Bertucci 2020). These KEGG pathways, such as glycolysis, which produces ATP from glycogen, and gluconeogenesis, which produces glucose from non-sugar organic matter, are all connected to muscle growth and development. In a refeeding experiment with rainbow trout, certain genes involved in glycolysis and gluconeogenesis were discovered (Johansen and Overturf 2006). Fatty acids are one of the primary energy sources for organisms since they may be oxidised and broken down to produce a substantial amount of energy (De Carvalho and Caramujo 2018). Proline and arginine are two amino acids that can be consumed to increase muscle growth and decrease extra fat (Wu 2013). In a certain range, increasing the level of faba beans in the feed increased the content of total saturated fatty acids, total monounsaturated fatty acids, polyunsaturated fatty acids and n-3 polyunsaturated fatty acids in the grass carp muscles. Exceeding the range would reduce its deposition in the muscles. The relationship between these KEGG pathways and muscle growth and development in fish can be further investigated in the future.

Studies have shown that material and energy metabolism can affect biological growth (Canosa and Bertucci 2020). This study also confirmed this point: the GO and KEGG enrichment results of DEGs mostly focused on material and energy metabolism. Six KEGG pathways involved in amino acid metabolism were shown to be up-regulated by feeding faba beans in the KEGG enrichment results, and there were alterations in the glycolysis/gluconeogenesis pathway. Fast-growing grass carp have been found to have dramatically increased expression of DEGs enriched in metabolic activities (Lu et al. 2020). These results indicated that the addition of faba beans had an important impact on the amino acid metabolism and energy metabolism of grass carp.

In this study, the construction of PPI network and the identification of hub genes were carried out for the DEGs enriched in different KEGG pathways linked to amino acid metabolism, fatty acid metabolism and immunity. A total of 10 hub genes were discovered, including ACOX3, ACADS, ACAT1, ECHS1, EHHADH, HADHA, HADH, HMGCL, MCCCC1 and MCCC2. During the feeding period of grass carp fed with faba beans, the immunity of crisp grass carp was reduced and its mortality increased. In this study, we found that many down-regulated genes in crisp grass carp were clustered into amino acid catabolism and immune-related category. For instance, down regulation of ECHS1 induces fatty acid and branched-chain amino acid BCAA accumulation (Qu et al. 2020). MCCC1 (Zhao et al. 2020) is a key enzyme involved in the BP of amino acid metabolism, which is related to the catabolism of valine and leucine. MCCC2 (Chen et al. 2021) also functions to regulate leucine metabolism.

ACOX3 has a role in regulating lipid homeostasis (Madureira et al. 2016). ACADS also regulate fatty acid metabolism and regulate energy homeostasis (Chen and Su 2015). ACAT1 (Abdelkreem et al. 2019) is involved in ketone body synthesis and degradation and isoleucine catabolism. Another study found that down regulation of ACAT1 was associated with muscle atrophy after spinal cord injury (Huang et al. 2021). EHHADH is also related to fatty acid metabolism. HADHA is a subunit of HADH that plays a role in fatty acid β-oxidation (Houten et al. 2012). HMGCL is the main enzyme for ketone body synthesis (Puisac et al. 2012).

The present study may indicate that the several hub genes screened from DEGs in the amino acid metabolism pathway are related to fatty acid metabolism, indicating that there is an interaction between fatty acid metabolism and amino acid metabolism to regulate the meat quality of grass carp, but the mechanism is still unclear and needs further research.

Conclusions

The outcomes of the current investigation revealed that dietary faba bean increased hardness, shear force and muscle fibre diameter in grass carp. These outcomes indicate that the DEGs of crisp grass carp were concentrated in material metabolism, energy metabolism and ribosome function, and the comprehensive regulation of multiple genes and multiple pathways changes the meat quality of crisp grass carp finally. The interaction between fatty acid metabolism and amino acid metabolism still needs further study.

Author contributions

Meilin Hao and Chuntao Li participated in the design of the study. Wenjie Cheng and Yuxiao Xie did the experiments and did the data analysis. Meilin Hao, Sumei Zhao and Wenjie Cheng drafted the manuscript and all authors contributed to finalising the writing.

Ethical approval

The experiment was carried out in accordance with the research plan of the Institutional Animal Care and Use Committee of Zunyi Normal College (protocol number: ZUNSHIFA[2018]08).

Consent form

All authors agree to participate in the publication of this article.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the fundings of this study are available from the corresponding author upon reasonable request.

Additional information

Funding

The research is supported by the Guizhou Province Colleges and Universities Youth Science and Technology Talent Development Project (Qian Jiao He KY[2020]106) and Zunyi Normal College PhD Start-Up Fund (BS[2019]26).