Abstract While tissue amino acid compositions reflect that of the dietary protein source, and the liver orchestrates amino acid metabolism. In this study, we investigated the muscle amino acid profiles in ordinary and crisp grass carp. The 22 amino acids were measured, and seventeen showed significant concentration differences. To understand the molecular mechanisms behind changes, we analyzed the liver transcriptome, and the 2519 differentially expressed genes (DEGs) were identified, with 1156 up-regulated and 1363 down-regulated genes. DEGs were enriched in ribosome-related biological processes. KEGG pathway analysis showed enrichment in tryptophan metabolism, lysine degradation, valine, leucine and isoleucine degradation, galactose metabolism, and glutathione metabolism with up-regulated genes, arginine and proline metabolism, arginine biosynthesis and alanine, aspartate, amino sugar and nucleotide sugar metabolism, N-Glycan biosynthesis and glutamate metabolism with down-regulated genes. A protein-protein interaction network with 260 nodes and 249 edges was constructed, and 3 modules were extracted. The top 10 hub genes with close connections to other nodes were ITM1, STT3B, SEL1L, UGGT1, MLEC, IL1B, ALG5, KRTCAP2, NFKB2, and IRAK3. In summary, this study identified candidate genes and focused on amino acid and glycan metabolism pathways, providing a reference for further investigation into liver amino acid metabolism in grass carp fed with broad beans. Keywords: Crisp grass carp, Broad bean, Amino acids, Transcriptomics 1. Introduction The quality of aquatic products is defined by a blend of their nutritional value and sensory characteristics [[33]1]. Both these aspects are subject to the influence of feeding patterns and dietary composition [[34]2,[35]3], which further enhance the overall characteristics of the products. Fish muscle is a highly sought-after protein source, boasting a wealth of vital amino acids. These amino acids not only serve as crucial building blocks for proteins but also play a pivotal role in the synthesis of other muscle components. Moreover, amino acids contribute significantly to the distinct flavor profile of the meat [[36]4]. The flavor constituents of the meat comprise both volatile and nonvolatile substances, as well as a variety of free amino acids. Amino acids are indispensable for the creation of diverse proteins that serve essential functions, such as transporting oxygen, producing vitamins, facilitating CO[2] exchange, catalyzing biochemical reactions, and forming structural proteins [[37]5]. Grass carp is one of the freshwater fish with the largest amount of cultivation in China, which also is the most consumed fish due to its delicious taste and rich protein content [[38]6]. Currently in fish feed, fishmeal and soybean meal are the main protein sources for grass carp [[39]7]. With the rapid development of the aquaculture industry, fishmeal and soybean meal became scarce and their prices increased [[40]8]. Therefore, it is essential to find an alternative plant protein feed for grass carp. The broad bean is a herbaceous plant belonging to the subfamily Papilionaceae of the family Fabaceae of the order Rosacea, and is widely planted all over the world. Broad beans are rich in protein, carbohydrates and trace elements [[41]9]. It is commonly used as a nutritional feed for animals such as pigs, poultry, ruminants and fish [[42]10]. After feeding a single feed of broad beans for 90–120 days, the muscle hardness increased, the meat was firm, and the taste was crisp [[43]6], [[44][11], [45][12], [46][13]], which is extremely popular among consumers. Compared to ordinary grass carp, crisped grass carp showed significant increase in muscle hardness, elasticity, chewing power and adhesion, collagen content, myofibril length and density, and reduction in myofibril diameter [[47]14,[48]15]. In addition, broad beans affect the fatty acid content of fish and increase fat deposition in viscera [[49]16]. It has been shown that continuous consumption of broad beans by grass carp leads to permanent inflammation-induced intestinal mucosal damage and hepatic steatosis [[50]17,[51]18], which seriously affects the health and quality of grass carp. The liver plays an important role in maintaining metabolic homeostasis as the main site of synthesis, metabolism and storage of carbohydrates, proteins and lipids [[52]19]. In order to investigate the potential mechanism by which amino acid composition is changed by fish feeding. This study quantified the muscle amino acids of ordinary and crisp grass carp, and pinpointed the amino acids exhibiting substantial variations among different groups. Additionally, we performed transcriptome sequencing analysis on the livers of grass carp and elucidated the impact of broad beans on the liver tissue metabolism of grass carp from a molecular biology standpoint. These findings serve as a valuable reference for promoting the well-being of grass carp breeding and enhancing their overall quality. 2. Materials and methods 2.1. Animal and sample collection 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 weighing approximately 768 ± 75 g at the beginning of the experiment, which were randomly distributed into two groups.: the crisp grass carp group (Group A) and the ordinary grass carp group (Group B), with three replicates for 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^−1; crude lipid: 43.8 g kg^−1; Tongwei Company, 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 the concentration of dissolved oxygen was ensured to exceed 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 and liver tissues, which was placed at −80 °C until RNA extraction. 2.2. Targeted quantitative analysis of amino acid The amino acid content of muscle tissue on ordinary and crisp grass carp was analyzed. The LC analysis was performed on ACQUITY Liquid chromatography (Waters, USA). Masss pectrometric detection of metabolites was performed on AB4000 (AB SCIEX, USA). L-Valine, L-Leucine, DL-Homocysteineand and L-Histidine were obtained from Sigma-Aldrich (Shanghai, China). L-Isoleucine were obtained from Aladdin (Shanghai, China). L-Aspar agine was obtained from OKA (Beijing, China). Glycine, L-Alanine, 4-Aminobutyric acid, L-Serine, L-Proline, L-Threonine, L-Ornithine hydrochloride, L-Aspartic acid, L-Lysine, L- Glutamine, Glutamic acid, L-Methionine, L-Phenylalanine, L-Arginine, L-Tyrosine and L-Tryptophan were all obtained from Sinopharm Chemical Reagent Co.Ltd. (Shanghai, China). DL-Tryptophan-2,3,3,-d3 was obtained from Pufen Biotechnology Co. LTD. (Pointe-Claire, Canada). Samples were extracted in 600 μL of 10 % formic acid in methanol-water(1:1, v/v). Add two steel balls and vortex for 30 s. Place in a tissue grinder and grind at 55 Hz for 90 s. Centrifuge at 12,000 rpm and 4 °C for 5 min. Take 20 μL of supernatant and add 580 μL 10 % formic acid in methanol-water(1:1, v/v)to dilute 30 times and vortex for 30 s. Take 100 μL of supernatant and add 100 μL Trp-d3(100 ng/mL), vortex for 30 s. The supernatant was filtered through 0.22 μm membrane, and the filtrate was added to the LC-MS bottle. ACQUITY UPLC® BEH C18 Column (2.1 × 100 mm, 1.7 μm, Waters, USA), The injection volume was 5 μL, the column temperature was 35 °C, and the mobile phase A-50 % methanol in water (containing 0.1 % formic acid), B-10 % methanol water (containing 0.1 % formic acid). The gradient elution protocol was as follows: from 0 to 6.5 min, 90-70 % B; from 6.5 to 7 min, 70-0% B; from 7 to 14 min, 0 % B; from 14 to 14.5 min, 0–90 % B; and from 14.5 to 17.5 min, 90 % B. The flow rate was set at 0.3 mL/min from 0 to 8.0 min, and at 0.4 mL/min from 8.0 to 17.5 min. The instrument used for ESI employed positive ionization mode. The ion source temperature was set at a scorching 500 °C, while the ion source voltage reached a whopping 5500 V. The collision gas was pumped up to 4 psi, while the curtain gas was cranked up to a remarkable 40 psi. To optimize the performance, both the atomizing gas and auxiliary gas were set to an intense 50 psi. The scans were meticulously conducted using the multiple reaction monitoring (MRM) technique, ensuring precise and reliable data collection. 2.3. RNA preparation Total RNA was isolated from liver tissue of ordinary and crisp grass carp using the TRIzol reagent (Takara, Dalian, China) according to the manufacturer's instruction. The concentration of the isolated RNA was determined by measuring absorbance at 260 nm. The integrity of the RNA was determined by agarose gel electrophoresis and Agilent BioAnalyzer 2100 (Agilent Technologies, San Jose, CA, USA). The RNA was used for transcriptomics analysis. 2.4. Library construction and sequencing Six RNA samples with high quality (RIN >8.7) and concentration (average concentration: 477.4 ng/μL) were used to construct the sequencing libraries and high-throughput sequencing was performed on the Illumina novaseq 6000 platform (San Diego, CA, USA) following the manufacturer's recommendations, generating 150 bp paired-end reads ([53]Table S1). The raw reads were effectively filtered and eliminated to obtain the high-quality clean reads: (1) the sequences containing adapters; (2) the sequences with more than 10 % of N bases; (3) the sequences with more than 50 % base quality values less than 10. 2.5. Differentially expressed genes analysis Illumina HiSeq 4000 sequencer reads were paired-end and quality controlled by Q30. Amplifification of the 30 adaptor and the removal of low-quality reads were performed by cutadapt software (v1.9.3), followed by alignment with the reference genome (C_idella_female_scaffolds. fasta V1) using hisat2 software (v2.0.4). Using the Ensembl gtf gene annotation file as a reference, the cuffdiff software was used to obtain the mRNA expression profiles and fold change, expressed as gene level fragments per kilobase per million (FPKM). The number of clean reads for each gene was calculated and FPKM was used to estimate the expression abundance of transcripts from different samples [[54]20]. Differential expression analyses of the A and B groups were performed using the DESeq R package [[55]21,[56]22] and genes with an p-value ≤0.05 and an expression | log2 Fold | ≥ 1were identified as DEGs. 2.6. GO and KEGG analysis To evaluate the function of these differentially expressed genes (DEGs), an analysis of Gene Ontology (GO) was performed using the GOseq software. The analysis covered three primary categories: biological processes, cellular components, and molecular functions. The KOBAS software was utilized to perform Kyoto encyclopedia of genes and genomes (KEGG) pathway analyses, identifying biological pathways that were enriched for the DEGs. 2.7. GSEA and PPI network analysis To further investigate the function of the identified genes, we performed Gene Set Enrichment Analysis (GSEA) to identify key pathways using KEGG annotation pathway sets. This analysis was conducted utilizing the clusterProfiler package and GSEABase (version 1.58.0) R software. Initially, the genes were ranked based on their expression levels. Subsequently, the enrichment scores were calculated according to this ranked-ordered gene list. Gene sets with a normalized enrichment score (NES) greater than 1 and a nominal P-value less than 0.05 were considered significantly enriched. To elucidate the interactive relationships among DEGs encoding proteins, we employed the STRING database ([57]https://cn.string-db.org/) for analysis and visualized the results using Cytoscape software (version 3.9.1). Subsequently, the highly correlated module was extracted from the entire PPI network using the Molecular Complex Detection (MCODE) algorithm via the Cytoscape MCODE plugin (version 2.0.2). Hub genes were identified using a degree centrality method in the Cytoscape plugin cytoNCA (version2.1.6), and the top 10 genes with the highest score were selected as hub genes. 2.8. qRT-PCR validation To validate the results of the transcriptome data, the top 10 hub genes were randomly selected, and their mRNA levels were measured using qRT-PCR. The specific primer sequences for qRT-PCR are as follows: Integral membrane protein 1 (ITM1, 5F: CGTGCCTGGATACATCTCCC; R: AATGGAGCCCGACTTCACTG), STT3 oligosaccharyltransferase complex catalytic subunit B (STT3B, 5F: GCGGAGATCGGCAATAAGGA; R: GGAGCTTGTGGTCCAGAGTC), SEL1L adaptor subunit of SYVN1 ubiquitin ligase (SEL1L, 5F: ACAGGGCTATGAGGTTGCAC; R: TTGATTCGGGCAACGGTGTA), UDP-glucose glycoprotein glucosyltransferase 1 (UGGT1, 5F: ATGGAGGGTCTGCGTAGTCT; R: TTGCACATTAGAGGGCCAGG), Malectin (MLEC, 5F: ATGGCTGTTGGTGGAAGTGT; R: GATAGGCAGACGCACACCAT), Interleukin-1 beta (IL1B, 5F: CCGCTCCACATCTCGTACTC; R: GTTTATGGAGCTGCCGGTCT), Asparagine-linked glycosylation 5 homolog (ALG5, 5F: CTGATCCTGGCTTTGGTCGT; R: CTGGGGAAATCCTCACGCTT), Keratinocyte-associated protein 2 (KRTCAP2, 5F: GGAGAACCTGGTCTTTGGCA; R: TTGTCACACAAACGCGATGC), Nuclear factor kappa B subunit 2 (NFKB2, 5F: GTACCTTGGAGGGAACTGGC; R: CATCCACAGGACCACCACTC), and Interleukin-1 receptor-associated kinase 3 (IRAK3, 5F: GGACCAACCTCACACTCTGG; R: GCATTCACATGGATCTGCCG). Total mRNA from the samples was extracted using the Trizol Kit. cDNA synthesis was then performed according to the manufacturer's instructions using the PrimeScript™ RT Reagent Kit with gDNA Eraser (Takara). Quantitative real-time PCR (qRT-PCR) was conducted using the Roche Step One real-time PCR system (Roche, USA). PCR amplification of all samples was performed in triplicate. β-Actin (5F: TGCCATGTATGTGGCCATCC; R: TCTTTCGGCTGTGGTGGTGA) was used as the normalization control to calculate the relative expression levels [[58]23]. 3. Results 3.1. Muscle amino acid profles in muscle of ordinary and crisp grass carp The quantitative analysis of 22 amino acids are shown in [59]Table 1. Among them, 13 amino acids (Phe, Tyr, Met, Val, Trp, Leu, Ile, Thr, Arg, Lys, Orn, Glu and Ser) were up-regulated, and 4 amino acids (Pro, His, Gln and Ala) were down-regulated in crisp grass carp compared with ordinary grass carp ([60]Table 2). ROC curve analysis was utilized to identify potential biomarkers. Various models were used to analyze the sensitivity values and area under the ROC curve (AUC) of the biomarkers that had been previously screened. Remarkably, 18 biomarkers exhibited an AUC ≥0.7 in all models. Among them, 17 amino acids showed high efectiveness (AUC >0.90) and 1 amino acids showed moderate efectiveness (AUC = 0.778) ([61]Fig. 1). As a result, we attained a relatively “good” level of accuracy by utilizing markers related to the metabolism of amino acids. Table 1. The quantitative analysis of 22 amino acids in muscle between ordinary and crisp grass carp. Amino acid(μg/g) Group_A1 Group_A2 Group_A3 Group_B1 Group_B2 Group_B3 Gly 690.4348178 811.5326733 749.060241 760.8379447 785.736 780.2857143 Ala 347.0064777 372.6463366 358.9951807 379.4158103 385.0992 379.8928571 GABA ND ND ND ND ND ND Ser 138.308502 149.7029703 140.2048193 85.24031621 88.0344 89.47857143 Pro 78.48582996 85.33782178 80.8626506 106.7976285 107.0136 107.0571429 Val 86.9611336 89.71485149 87.18072289 42.83715415 43.74072 44.23 Thr 162.5611336 163.8392079 159.860241 128.3620553 128.0232 127.0714286 Ile 82.82186235 88.69544554 85.76385542 36.23833992 36.29376 36.24214286 Leu 134.562753 137.0352475 136.6481928 69.71762846 71.17704 74.85 Asn 14.32639676 15.44435644 13.28819277 13.81588933 17.35416 12.94785714 Orn 103.9846154 110.2170297 109.7421687 61.80474308 64.494 66.635 Asp 11.22267206 12.5230099 11.26409639 11.5484585 12.80592 11.81714286 Hcy ND ND ND ND ND ND Gln 115.7902834 117.8732673 120.0144578 141.5525692 150.1776 147.7357143 Lys 623.9951417 658.5077228 663.2385542 358.3351779 360.5904 364.8357143 Glu 26.14736842 28.02011881 26.54819277 15.77454545 15.71976 16.55142857 Met 51.54048583 52.66336634 52.13783133 30.67897233 30.5172 31.55214286 His 433.4866397 454.7762376 457.4891566 1008 1074.096 1089.214286 Phe 60.0048583 59.11485149 60.11710843 25.76703557 26.64936 27.26285714 Arg 208.8510121 221.2681188 221.4939759 96.28221344 100.5192 101.5785714 Tyr 75.88421053 77.60316832 78.03614458 42.63438735 42.44184 43.13357143 Trp 22.70987854 23.29663366 23.72674699 9.556363636 9.63288 10.07357143 [62]Open in a new tab Note: ND indicates that the device is not detected. Table 2. List of statistically signifcant diferences in amino acid concentrations between ordinary and crisp grass carp (n = 3). Amino acid Mean_A SD_A Mean_B SD_B Fold Change P value Phe 59.74563333 0.549104374 26.55976667 0.751967289 2.2495 4.12E-07 Tyr 77.1745 1.138203044 42.7366 0.357044031 1.8058 9.57E-07 Met 52.1139 0.561831389 30.9161 0.556701724 1.6857 1.29E-06 Val 87.95223333 1.530457897 43.60263333 0.706590237 2.0171 1.39E-06 Trp 23.2444 0.510405907 9.7543 0.27915485 2.383 2.30E-06 Leu 136.0820667 1.32987618 71.91486667 2.64456383 1.8923 3.00E-06 Ile 85.7604 2.936751564 36.25806667 0.031004247 2.3653 8.20E-06 Thr 162.0868333 2.031454381 127.8189 0.669164023 1.2681 1.00E-05 Arg 217.2043667 7.235109447 99.46 2.802568772 2.1838 1.25E-05 His 448.584 13.14491184 1057.103433 43.19144733 0.42435 2.00E-05 Lys 648.5804667 21.42254757 361.2537667 3.300631237 1.7954 2.13E-05 Orn 107.9812667 3.469346782 64.31123333 2.420331024 1.679 5.75E-05 Glu 26.90523333 0.986081094 16.01523333 0.465138736 1.68 6.55E-05 Ser 142.7387667 6.105276232 87.58443333 2.154680771 1.6297 0.000122797 Pro 81.5621 3.479130102 106.9561 0.138977516 0.76258 0.000226111 Gln 117.8927 2.112166821 146.4886333 4.445676293 0.80479 0.000548486 Ala 359.5493333 12.8288789 381.4693 3.152623766 0.94254 0.045288128 [63]Open in a new tab Note: Mean_∗: average value; SD_∗: standard deviation. Fig. 1. [64]Fig. 1 [65]Open in a new tab The top list of amino acids selected as potential biomarkers by ROC curve analysis. AUC (0.5–0.7), low accuracy; AUC (0.7–0.9), moderate accuracy; AUC (>0.9), high accuracy. In order to gain a deeper insight into the metabolic distinctions between ordinary and crisp grass carp for individual amino acids, we used cluster heatmaps to visualize the identified distinct metabolites. The majority of samples exhibited a distinct separation into two distinct clusters ([66]Fig. 2). Fig. 2. [67]Fig. 2 [68]Open in a new tab Heat maps for identified metabolites in ordinary and crisp grass carp. The color of each section is proportional to the signifcance of change of metabolites (red, upregulated; blue, downregulated). Rows, samples; columns, amino acid. 3.2. Differentially expressed genes of liver tissue between ordinary and crisp grass carp In this study, we used high-throughput RNA sequencing to construct six libraries in liver tissue. The number of clean reads for each sample varied between 35 million and 43 million, demonstrating the robustness of our data. Moreover, we achieved an impressive mapping rate of over 92.13 %, with a remarkable 94.16 % of reads mapping uniquely. In addition, more than 84.12 % mapped to gene were detected in each sample ([69]Table 3). Based on the correlation analysis performed on 6 samples, it was observed that the variations among biological replicates were minimal and the consistency was remarkably high. This finding strongly suggests that the experimental samples was both dependable and uniform ([70]Fig. 3A and B). Upon comparing the liver transcriptomes of ordinary and crisp grass carp, a total of 2519 differentially expressed genes (DEGs) were identified. In order to meet the selection criteria, these DEGs had to satisfy the conditions of |log2FoldChange| > 1 and a P-value of less than 0.05. Among these DEGs, 1156 were found to be up-regulated, while 1363 were down-regulated ([71]Fig. 3C). However, after adjusting the P-value using the p-adjust method with a threshold of less than 0.05, 1399 genes were detected in liver tissues. It is worth noting that a significant number of these differentially expressed genes exhibited high expression levels in the liver tissues of two groups. Table 3. Summary of reads and matches. Samples Clean Reads Mapped Reads Unique Mapped Reads Mapped to Gene Group_A1 36354078 92.70 % 94.27 % 84.78 % Group_A2 39272704 92.13 % 94.34 % 84.12 % Group_A3 35914404 92.85 % 94.16 % 85.50 % Group_B1 42460834 93.83 % 94.73 % 86.63 % Group_B2 40967930 93.82 % 94.80 % 86.47 % Group_B3 39474196 93.81 % 95.01 % 86.52 % [72]Open in a new tab Fig. 3. [73]Fig. 3 [74]Open in a new tab Transcriptome profiles and DEGs analysis. (A) The correlation analysis of 6 samples. (B) Principal components analysis (PCA) scores plot. (C) Volcano plot of differentially expressed genes (DEGs) (D–F) The significant enrichment of GO pathway for DEGs. (G–I) The significant enrichment of KEGG pathway for DEGs. To dissect the functional categories of differentially expressed genes (DEGs), Gene Ontology (GO) enrichment analysis was performed. According to the GO enrichment analysis results, the majority of the differentially expressed genes (DEGs) were categorized into three primary functional groups, namely cellular composition, biological progression, and molecular capability ([75]Fig. 3D). In the cellular component category, most genes were enriched in nuclear outer membrane-endoplasmic reticulum membrane network (GO:0042175, 53 genes), endoplasmic reticulum membrane (GO:0005789, 52 genes), endomembrane system (GO:0012505, 118 genes), endoplasmic reticulum (GO:0005783, 40 genes) and organelle membrane (GO:0031090, 96 genes). In the molecular function category, DEGs were enriched in iron ion binding (GO:0005506, 28 genes), oxidoreductase activity (GO:0016491, 72 genes), heme binding (GO:0020037, 24 genes), tetrapyrrole binding (GO:0046906, 24 genes) and catalytic activity (GO:0003824, 339 genes). Meanwhile, the DEGs involved in fatty acid metabolic process (GO:0006631, 21 genes), small molecule metabolic process (GO:0044281, 78 genes), lipid metabolic process (GO:0006629, 54 genes), and oxidation-reduction process (GO:0055114, 69 genes) were enriched in the biological process category ([76]Table S2). GO enrichment showed that up-regulated genes were significantly enriched for cellular components and biological process ([77]Fig. 3E). In the molecular function category, up-regulated genes were enriched in iron ion binding (GO:0005506, 19 genes), oxidoreductase activity (GO:0016491, 42 genes), monooxygenase activity (GO:0004497, 16 genes), heme binding (GO:0020037, 16 genes), tetrapyrrole binding (GO:0046906, 16 genes), and oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen (GO:0016705, 17 genes). Meanwhile, in the biological process category, up-regulated genes were involved in oxidation-reduction process (GO:0055114, 38 genes) mainly ([78]Table S3). GO enrichment showed that down-regulated genes were significantly enriched for cellular components ([79]Fig. 3F). In the molecular function category, down-regulated genes were enriched in endomembrane system (GO:0012505, 95 genes), nuclear outer membrane-endoplasmic reticulum membrane network (GO:0042175, 45 genes), endoplasmic reticulum membrane (GO:0005789, 44 genes), endoplasmic reticulum (GO:0005783, 37 genes), and organelle membrane (GO:0031090, 74 genes) ([80]Table S4). To identify the primary biochemical metabolic pathways and signal transduction pathways, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway classification and functional enrichment analysis were conducted for the differentially expressed genes (DEGs). According to the KEGG pathway enrichment analysis, it was observed that the differentially expressed genes (DEGs) exhibited significant enrichment in a total of 20 pathways, as depicted in [81]Fig. 3G. Notably, 10 of these pathways were found to be closely associated with various aspects of metabolism, such as energy metabolism, amino acid metabolism, carbohydrate metabolism, and lipid acid metabolism. Specifically, these pathways included steroid biosynthesis, glutathione metabolism, metabolism of xenobiotics by cytochrome P450, fatty acid degradation, drug metabolism - cytochrome P450, terpenoid backbone biosynthesis, N-Glycan biosynthesis, and steroid hormone biosynthesis. The genes in these pathways were up-regulated ([82]Table S5). According to the KEGG pathway enrichment analysis, it was found that the genes were up-regulated exhibited significant enrichment in 20 pathways ([83]Fig. 3H). Notably, 17 of these pathways were closely associated with metabolism and organismal systems. The metabolism pathways encompassed various processes such as xenobiotic metabolism by cytochrome P450, drug metabolism-cytochrome P450, fatty acids degradation, arachidonic acid metabolism, drug metabolism-other enzymes, linoleic acid metabolism, tryptophan metabolism. On the other hand, the organismal systems pathways included protein digestion and absorption, longevity regulating pathway-worm, cholesterol metabolism, ovarian steroidogenesis, vitamin digestion and absorption, and carbohydrate digestion and absorption of carbohydrates ([84]Table S6). The KEGG pathway enrichment analysis revealed a significant enrichment of down-regulated genes in 20 pathways ([85]Fig. 3I). Notably, 13 of these pathways were found to be closely associated with metabolism and organismal systems. Metabolism pathways included amino sugar and nucleotide sugar metabolism, steroid biosynthesis, terpenoid backbone biosynthesis, N-Glycan biosynthesis, various types of N-glycan biosynthesis, fatty acid elongation, and fatty acid biosynthesis. Organismal systems pathways included IL-17 signaling pathway, fat on and absorption, adipocytokine signaling pathway, osteoclast differentiation, and C-type lectin receptor signaling pathway ([86]Table S7). This result, in consistence with GO analysis, further indicating metabolic and immune enhancement in crisp grass carp. 3.3. Delineation of key genes via GSEA and PPI network analysis To globally assess the diferences between the ordinary and crisp grass carp liver transcriptome, we conducted Gene Set Enrichment Analysis (GSEA) to identify enriched gene sets, emphasizing both positive and negative regulation as per KEGG functional annotations. Detailed results can be found in [87]Tables S8 and S9. A total of 32,818 genes ranked by fold change in liver gene expression were subjected to Gene Set Enrichment Analysis (GSEA). The Enrichment Score (ES) reflects the degree to which a gene set is overrepresented at the top or bottom of the ranked list of genes. A positive ES indicates upregulation, while a negative ES indicates downregulation of a specific gene set. This study focused metabolic pathways, which could deepen our understanding of the mechanism of meat quality regulation. GSEA for KEGG pathway mapping showed that these pathways were mainly involved in amino acid metabolism, lipid metabolism, Carbohydrate metabolism and aenobiotics biodegradation and metabolism, the most upregulated pathways in the ordinary grass carp liver were terpenoid backbone biosynthesis (ES = −0.75), amino sugar and nucleotide sugar metabolism (ES = −0.62), N-Glycan biosynthesis (ES = −0.51),Various types of N-glycan biosynthesis (ES = −0.43) fatty acid biosynthesis (ES = −0.45), fat digestion and absorption (ES = −0.35) and steroid biosynthesis (ES = -0.71) ([88]Fig. 4A), whereas the most downregulated gene sets belonged to folate biosynthesis (ES = 0.63), tryptophan metabolism (ES = 0.68), lysine degradation (ES = 0.60), fatty acid degradation (ES = 0.63), valine, leucine and isoleucine degradation (ES = 0.61) and drug metabolism - cytochrome P450 (ES = 0.51) ([89]Fig. 4B). Fig. 4. [90]Fig. 4 [91]Open in a new tab GSEA functional analysis between the ordinary and crisp grass carp. The enrichment plots using KEGG as the database showed the most enriched KEGG terms of the upregulated (A) or downregulated (B) gene sets in the ordinary grass carp. The metabolic pathway of differential gene enrichment in crisp grass carp liver tissues was further analyzed. Matebolism pathway enrichment analysis of the up-regulated genes including Amino acid metabolism, Lipid metabolism, Metabolism of cofactors and vitamins and xenobiotics biodegradation and metabolism ([92]Fig. 5A). Matebolism pathway enrichment analysis of the down-regulated genes including amino acid metabolism, carbohydrate metabolism, and lipid metabolism ([93]Fig. 5B). After conducting STRING analysis on the differentially expressed genes (DEGs), we proceeded to build and visualize the protein-protein interaction (PPI) network using Cytoscape. The resulting network consisted of 260 nodes and 249 interactions, as shown in [94]Fig. 6A. To further analyze the entire PPI network, we utilized cytoHubba. After evaluating the connectivity level of every individual node, the top 10 hub genes with the strongest associations to other nodes were identified as ITML, STT3B, SEL1L, UGGT1, MLEC, IL1B, ALG5, KRTCAP2, NFKB2 and IRAK3. The whole PPI network was analyzed by MCODE ([95]Fig. 6B). We extracted a total of 8 modules from the PPI network. Subsequently, for functional analysis, we specifically identified three modules (Modules 1–3) that exhibited a significant MCODE score (>3) as well as a substantial number of nodes (>3). The pathway enrichment analysis uncovered that the differentially expressed genes (DEGs) in module 1 exhibited predominant enrichment in glycan biosynthesis and metabolism. The hub genes ITM1, STT3B, SEL1L, UGGT1 and MLEC participated in the pathway. Module 2 was mainly associated with immunity, and included the hub genes IL1B, NFKB2, BCL3, and IRAK3. Module 3, including RAB7, STX8, RAB20, and TBC1D2, exhibited a close relationship with intrinsic transport and secretion processes of the cell. After STRING analysis of the up-regulated genes, the PPI network was constructed and analyzed by cytoHubba ([96]Fig. 6C). After the connectivity degree of each node was calculated, BBOX1 were the hub gene with the closest connections to other nodes. After STRING analysis of the down-regulated genes, the PPI network was constructed and analyzed by cytoHubba ([97]Fig. 6D). After the connectivity degree of each node was calculated, NANSA, ALG5, NSDH1 were the hub genes with the closest connections to other nodes. Fig. 5. [98]Fig. 5 [99]Open in a new tab The most significant KEGG signaling pathway related with matebolism in crisp grass carp. (A) upregulated pathway. (B) downregulated pathway. Fig. 6. [100]Fig. 6 [101]Open in a new tab Identification of key metabolic pathways and hub genes. (A) The PPI network constructed and visualized by Cytoscape, including 260 nodes and 249 edges. (B) Module analysis. Blue nodes represent the DEGs, red nodes represent the genes upregulated and the green nodes represent the genes upregulated in crisp grass carp, the lines represent the interaction between two nodes. (C) Matebolism pathway enrichment analysis of the up-regulated genes. (D) Matebolism pathway enrichment analysis of the down-regulated genes. (E) Validation of the expression of 10 genes in transcriptome data by RT-qPCR. Relative mRNA expression level was calculated using 2^−^ΔΔCt method. All experiments were performed using three biological replicates. Group A: crisp grass carp, Group B: ordinary grass carp. To validate the RNA-Seq results obtained, 10 hub genes were randomly selected for validation by qRT-PCR. The qPCR expression patterns of the selected genes, as illustrated in [102]Fig. 6E, were consistent with the RNA-Seq analysis results. This alignment underscores the reliability and accuracy of the RNA-Seq methodology employed in this study. 4. Discussion Aquaculture provides people with high-quality protein [[103]24]. The crisp grass carp is very popular among consumers because of its improved taste, and some products are exported to Southeast Asia and North America [[104]13], [[105]25], [[106]26]. The meat sensory characteristics and taste are enhanced by the degradation of peptides and amino acids, such as anserine and carnosine, thereby leading to an improved overall flavor profile [[107]4]. The maillard reaction is responsible for the creation of most flavor substances found in food. This intricate chemical process occurs when carbonyl and amino compounds, like reducing sugars and amino acids, undergo a series of reactions, resulting in the formation of diverse and aromatic flavor compounds. The delicious taste associated with meat is contributed by aspartic acid, glutamic acid, and various salts of glutamine, asparagine, taste peptides, and nucleotides, lending it a distinct and refreshing palate. In the present study, 17 amino acid differences between ordinary and crisp grass carp. Among them, 13 amino acids were up-regulated, and 4 amino acids were down-regulated in crisp grass carp. The result indicate that amino acids play a crucial role in providing essential nutrients and enhancing the flavor of crisp grass carp. Glutamine and alanine play crucial roles in transporting carbon and nitrogen between tissues as part of the glucose-alanine and glucose-glutamine cycles, and they serve as key substrates for gluconeogenesis. Arginine, an indispensable amino acid for protein synthesis, has been found to potentially enhance intramuscular fat (IMF) synthesis by increasing the expression of vital lipogenesis genes in muscle tissues [[108]27]. Additionally, the presence of low levels of aromatic amino acids, particularly proline, has been associated with significantly stronger antioxidant properties, as demonstrated by their potent 2,2-diphenyl-1-picrylhydrazyl (DPPH) radical scavenging activity [[109]28]. In the present study, databases of transcriptomes were established by utilizing liver samples from grass carp. One reason was that the liver plays an important role in maintaining the metabolic stability of fish [[110]19]. Another reason was that dietary nutrient levels or dietary changes have significant effects on the liver of fish. Studies have shown that broad bean affects the amino acids and fatty acids content of grass carp and further affects the meat quality of fish [[111]6]. The number of clean reads in each sample varied between 35 million and 43 million, more than 84.12 % mapped to gene were detected in each sample. Our results also showed that the sequencing depth is sufficient and the sequencing quality is high enough to meet the requirements of later analysis. The 2519 DEGs that have been identified provide valuable insights into the contrasting metabolic profiles between ordinary and crisp grass carp. These distinct metabolic differences highlight the importance of further research on the liver metabolism mechanism in crisp grass carp. In the present study, GO and KEGG analysis results of DEGs showed that the liver metabolism capacity of crisp grass carp were enhanced. Among them, lipid acid metabolism, amino acid metabolism, glycan biosynthesis and metabolism were significantly enriched. This may also further explain the characteristics of crisp grass carp. The liver mainly performs metabolic functions in the body [[112]29]. Previous studies have shown that crisp grass carp has hard meat and crisp taste, which may be caused by the enhancement of liver lipid acid metabolism and amino acid metabolism [[113]30,[114]31]. Studies have shown that feeding broad beans enhances muscle stiffness and liver metabolism in fish, which is consistent with the results of this study [[115]32,[116]33]. Using GSEA, we have demonstrated that the amino sugar and nucleotide sugar metabolism, and the tryptophan metabolism, lysine degradation and valine, leucine and isoleucine degradation were up-regulation, whereas the amino sugar and nucleotide sugar metabolism, steroid biosynthesis, N-Glycan biosynthesis, Various types of N-glycan biosynthesis were down-regulation in crisp compared with ordinary grass carp liver tissue, it seems that there was a metabolic reprogramming toward increased amino acid metabolism but reduced amino sugar and nucleotide sugar metabolism in crisp grass carp. Moreover, the hub genes ITM1, STT3B, SEL1L, UGGT1 and MLEC are mainly involved in the regulation of glycan biosynthesis and metabolism. These glycosyl biosynthetic and metabolic pathways play critical roles in the normal function of cells and organisms [[117]34]. For example, the synthesis and modification of glycoproteins and glycolipids are critical for cellular signaling, cell adhesion, and cell-cell interactions. UGGT1 is an enzyme in the endoplasmic reticulum (ER), which is involved in the glycosylation modification process of proteins [[118]35]. ALG5 is a glycosyltransferase involved in the synthesis of N-glycan chains of proteins. ALG5 and UGGT1 genes play a crucial role in the protein quality control of the endoplasmic reticulum and the maintenance of cellular homeostasis. IL1B (Interleukin-1 beta) is a cytokine (cytokine), which plays an important regulatory role in the immune system and inflammatory response. In this study, IL1B and IRAK3 genes down-regulated in liver tissue of the crisp grass carp. A series of inflammatory reactions in the body will lead to excessive levels of IL1B [[119]36]. According to recent research conducted by Rowley et al. [[120]37], the role of Interleukin-1 receptor-associated kinase 3 (IRAK3) as a pseudokinase mediator in the human inflammatory pathway has been extensively studied. Interestingly, disrupting its function has been found to significantly boost the body's natural ability to combat tumor growth, indicating a potential link between IRAK3 ablation and improved antitumor immune response. The study observed a reduction in IL1B transcripts of ordinary campared with crisp grass carp, suggesting that Broad beans also contain undesirable factors such as trypsin inhibitors, which can lead to an increase in free radicals and subsequently cause oxidative stress. These factors contribute to the development of an inflammatory environment. KRTCAP2 (Keratinocyte associated protein 2) is an intracellular structural protein that is related to keratin and may play a role in cellular structures such as the cytoskeleton and intercellular connections [[121]38,[122]39]. NFKB (nuclear factor kappa B) is the important gene in the nuclear factor kappa B pathway. The pathway of nuclear factor kappa B serves as a crucial cell signaling pathway, exerting significant influence on the immune and inflammatory responses of cells [[123][40], [124][41], [125][42]]. These hub genes were significantly up-regulated in the liver of crisp grass carp, confirming the results of the GO and KEGG pathways. Feeding broad beans enhanced the metabolic capacity of grass carp liver, especially glycan biosynthesis and metabolism. Accordingly, aspartic acid, glycine, and glutamic acid are known to promote the healing of wounds [[126]43]. On the other hand, tyrosine, methionine, histidine, lysine, and tryptophan exhibit remarkable efficacy in neutralizing radicals during oxidative reactions [[127]44]. Additionally, hydrophobic amino acids have the ability to interact with membrane lipid bilayers to effectively target and eliminate radicals. Therefore, in this study, several of these potential genes are linked to the immune system, indicating that the intake of fish muscle protein is correlated with advantageous health effects, notably in terms of reducing inflammation, enhancing antioxidant activity, and combating microbial infections. The results were highly consistent with other studies on fish [[128]45,[129]46], and more changes in their functions and immune mechanisms still need to be further investigated. In the previous study, we have performed transcriptome sequencing on the muscle tissue of ordinary and crisp grass carp and to screen out differentially expressed genes and signaling pathways [[130]47]. It was shown that amino acid metabolism were closely related to meat quality, 6 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 [[131]47]. Moreover, several hub genes screened from DEGs in the amino acid metabolism pathway are related to fatty acid metabolism and Glycolysis/Gluconeogenesis [[132]47], indicating that there is an interaction among fatty acid metabolism, carbohydrate and amino acid metabolism to regulate the meat quality of grass carp, but the mechanism is still unclear. In the present study, we also found that these amino acid metabolism-related pathways play an important role in liver tissue, and the identified candidate genes focused on glycan biosynthesis and metabolism pathways. This seems to effectively explain part of the results of previous studies. Amino acids are widely involved in metabolism, are the basic building blocks of peptides proteins and substrates for gluconeogenesis. Hepatic uptake is remarkably efficient for most amino acids, with the notable exception of branched-chain amino acids. This efficiency underscores the liver's pivotal role in orchestrating amino acid metabolism and regulating systemic amino acid levels. 5. Conclusions Our study first revealed that the incorporation of broad bean induced significant changes in the hepatic transcriptome and amino acid profiles of muscle in grass carp. There was a metabolic reprogramming toward increased amino acid metabolism but reduced amino sugar and nucleotide sugar metabolism in crisp grass carp, PPI analysis revealed that the network modules were tightly correlated with glycan biosynthesis and metabolism changes, several hub genes of which were proposed to be crucial for amino acid metabolism, which further indicates that broad bean diet enhances liver metabolism and immune response, especially glycan biosynthesis and metabolism. This study might provide some insights into the underlying liver metabolic changes in fish fed with a broad bean diet. Data availability statement The datasets analyzed during the current study are available in the [NCBI BioProject repository, PRJNA1053764]. The preprint version to this article can be found online at [133]https://doi.org/10.21203/rs.3.rs-3320206/v1. Funding This work was supported by Guizhou Province colleges and universities youth science and technology talent development project (Grant numbers [Qian Jiao He KY[2020]106]) and Zunyi Normal College PhD start-up fund (Grant numbers [BS[2019]26]). Ethics statement All experimental procedures in this study strictly adhered to the regulations outlined by Chinese law regarding animal research. Furthermore, the comprehensive protocol received official approval from the Ethics Scientific Committee for Experiments on Animals at Zunyi Normal College (Protocol number: ZUNSHIFA[2018]08). Abbreviations GO Gene Ontology ALG5 Asparagine-linked glycosylation 5 homolog DEGs Differentially expressed genes ER Endoplasmic reticulum GSEA Gene Set Enrichment Analysis IL1B Interleukin-1 beta IRAK3 Interleukin-1 receptor-associated kinase 3 ITM1 Integral membrane protein 1 KEGG Kyoto encyclopedia of genes and genomes KRTCAP2 Keratinocyte associated protein 2 MLEC Malectin NFKB Nuclear factor kappa B PPI Protein-protein interaction SEL1L SEL1L adaptor subunit of SYVN1 ubiquitin ligase STT3B STT3 oligosaccharyltransferase complex catalytic subunit B UGGT1 UDP-glucose: glycoprotein glucosyltransferase 1 [134]Open in a new tab CRediT authorship contribution statement Meilin Hao: Writing – review & editing, Validation, Supervision, Software, Project administration, Investigation, Funding acquisition, Formal analysis. Junhong Zhu: Writing – original draft, Formal analysis. Yuxiao Xie: Software, Formal analysis. Wenjie Cheng: Formal analysis. Lanlan Yi: Formal analysis. Sumei Zhao: Writing – review & editing, Visualization, Resources, Project administration, Data curation. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Footnotes ^Appendix A Supplementary data to this article can be found online at [135]https://doi.org/10.1016/j.heliyon.2024.e38323. Appendix A. Supplementary data The following are the Supplementary data to this article: Multimedia component 1 [136]mmc1.xlsx^ (18.8KB, xlsx) Multimedia component 2 [137]mmc2.xlsx^ (1.8MB, xlsx) References