Transcriptomic assessment of dietary fishmeal partial replacement by soybean meal and prebiotics inclusion in the liver of juvenile Pacific yellowtail (Seriola lalandi)

Seriola lalandi is an important species for aquaculture, due to its rapid growth, adaptation to captivity and formulated diets, and high commercial value. Due to the rise in fishmeal (FM) price, efforts have been and still are made to replace it partially or entirely with vegetable meals in diets for carnivorous fish. The use of prebiotics when feeding vegetable meals has improved fish health. Four experimental diets were assessed in juveniles, the control diet consisted of FM as the main protein source, the second diet included 2% of GroBiotic®-A (FM-P), in the third diet FM was partially replaced (25%) by soybean meal (SM25), and the fourth consisted of SM25 with 2% of GroBiotic®-A (SM25-P). Growth was evaluated and RNA-seq of the liver tissue was performed, including differential expression analysis and functional annotation to identify genes affected by the diets. Growth was not affected by this level of FM replacement, but it was improved by prebiotics. Annotation was achieved for 59,027 transcripts. Gene expression was affected by the factors: 225 transcripts due to FM replacement, 242 due to prebiotics inclusion, and 62 due to the interaction of factors. The SM25-P diet showed the least amount of differentially expressed genes against the control diet. The replacement of FM (25%) by soybean meal combined with prebiotics (2%) represents a good cost-benefit balance for S. lalandi juveniles since the fish growth increased and important metabolic and immune system genes in the liver were upregulated with this diet.


Introduction
The Pacific yellowtail Seriola lalandi is a pelagic marine fish of the perciform order of the Carangidae family, which is distributed worldwide in temperate and subtropical waters [1] ranging from 18 to 24 °C [2,3]. It feeds mainly on mackerel, anchovy, sardine, and squid [2] and can grow up to 2.5 m long and weigh up to 70 kg [4]. It is an important species for commercial and recreational fishing in Japan, Australia, and New Zealand [5]. S. lalandi is a commercially important species for aquaculture because it grows rapidly, adapts well to captivity, accepts formulated diets, and has an excellent acceptance in the market and high commercial value due to the quality of its meat [6]. One of the key aspects for the successful culture of this species is adequate feeding and nutrition, which depends on the quality and source of the ingredients that make up the diet [2]. Fishmeal is the primary source of protein in diets for carnivorous fish. This ingredient has a high digestibility and an excellent amino acid profile [7]. However, because it is a limited resource, in the last years the costs have increased dramatically, and its availability has decreased since its global production has not increased to meet the demands of the aquaculture industry [8]. Replacing fishmeal with alternative protein sources helps to optimize production costs, but also it may reduce the fishing pressure on the species that are used to produce the fishmeal, thus reducing the impact on marine ecosystems [9]. In this sense, several studies have been performed to replace fishmeal with vegetable meals [10,11]. Soybean meal (SBM) has a good balance of essential amino acids and an appropriate nutritional composition at a reasonable price, it is widely used to replace fishmeal [12,13]. However, it should be considered that SBM has several antinutritional factors (ANFs) [14]. Besides, when fed to carnivorous fish SBM can cause intestinal damage leading to enteritis [10] and the ANFs produced by plants as a defense mechanism against animals can alter the biochemical, physiological, or immunological response of the animals fed with SBM [15][16][17].
The use of prebiotics in the diets of fish has emerged as an alternative to reduce the adverse effects of alternative proteins such as SBM-based formulated diets [18][19][20]. Some commercial products containing prebiotics have recently been evaluated with interesting results in fish such as Pre-vida™, Bio-MOS®, and GroBiotic®-A [21][22][23]. Prebiotics are indigestible carbohydrates used as food ingredients that beneficially affect the host, by stimulating the growth and activating the metabolism of the health-promoting bacteria in the digestive tract [24]. GroBiotic®-A is a widely used prebiotic supplement in aquaculture, which consists of a mixture of partially autolyzed brewer's yeast and dry fermentation products. It is composed of several polysaccharides, one of which, β-glucan is known to induce positive immunological responses in fish [25]. It has been reported that in juvenile hybrid striped bass (Morone chrysops × Morone saxatilis) the inclusion of GroBiotic®-A in the diet, improved respiratory burst and resistance against infection with Streptococcus iniae and showed a significant increase in feed efficiency [19]. In another work, [20] authors reported an improved growth performance in hybrid striped bass fed a diet enriched with GroBiotic®-A and a general improvement in survival. Likewise, it was reported that when GroBiotic®-A was included in the diet, the apparent digestibility coefficients of protein and carbohydrate increased for red drum (Sciaenops occelatus) [18]. Other works have shown that prebiotics affect lipid and glucose metabolism in mammals, reducing hepatic lipogenesis, triglycerides levels, and cholesterol in serum and liver while improving glucose tolerance [26,27].
Differences in diet composition promote changes in the levels of enzyme activities involved in the intermediary metabolism in the fish liver [28]. Metabolic factors, such as metabolite and enzymatic activity, are currently used to determine the capacity for metabolic adaptation to dietary supply in fish. In this regard, key liver enzyme activities of intermediary metabolism have been shown to match well with nutritional status and the growth rate of fish [29]. The liver is the primary metabolic organ of fish, being the center of the intermediary metabolism [29,30]. It is involved in the use of nutrients, plays an essential role in carbohydrate and lipid metabolism, synthesis of bile salts and most of the plasma proteins and hormones, and is the main energy storage. These features make the liver a reliable indicator in the evaluation of the nutritional condition, including energy metabolism, and an indicator of the state of the fish immune system [31][32][33][34]. Moreover, the liver has been the target organ to assess the effect of experimental diets on the physiological performance of fish [29,30,32,[35][36][37][38][39][40][41][42][43][44]; and recently, there is a growing body of research on liver transcriptomic responses to alternative diets for fish [30,[45][46][47][48][49].
In this study, the effects of three experimental diets on the liver physiology of juveniles were evaluated using a transcriptomic approach by implementing the RNA-seq technology. The present work aims to identify which genes were affected in their expression levels due to the FM replacement by SBM, the prebiotic inclusion, and their combination, making emphasis on key biological processes like carbohydrate metabolism, lipid metabolism, growth, and the immune response.

Experiment design and feeding trial
The juveniles of Pacific yellowtail S. lalandi used in this work came from Ejido Eréndira, Ensenada, B.C., Mexico, which were donated from the company Ocean Baja Labs. The feeding trial was conducted in the aquaculture department of the Center for Scientific Research and Higher Education of Ensenada (CICESE). Before the experiment, the fish were placed for 17 days in 3 m 3 raceway ponds, then reared for 58 days in a 9 m 3 pond where they reached the size of 97.9 ± 8.5 g. For the experimental trial, 180 fish were randomly taken and distributed in 12 tanks at a rate of 15 organisms per tank. The tanks consisted of circular fiberglass ponds with a capacity of 500 L, with a flat bottom, a biofilter with Kapdes type aeration, with a volume of 350 L, a compensation tank of 275 L, and two water pumps (Sweetwater High efficiency She 3.0 Aquatic Eco-Systems and an Aqua Logic heat pump model 2TWB0018A1000AB). In the experimental tanks, dissolved oxygen was maintained between 5 and 7 mg/L, the temperature at 23.5 ± 1 °C, and total ammoniacal nitrogen below 0.05 mg/L. The fish were in a 12-day acclimatization period in the new culture system, where they fed the commercial diet for marine fish (Skretting: 55% protein and 15% lipids), three times a day (8:00, 12:00, and 16:00 h). The feeding trial was conducted to evaluate the effect of including a commercial prebiotic, GroBiotic®-A, and the replacement of fishmeal with SBM. Four experimental diets were formulated (Table 1), the control diet consisted of fishmeal as the main protein source (FM), the second diet was formulated with fishmeal as the main protein source with the inclusion of 2% GroBiotic®-A (FM-P), the third diet included fishmeal with 25% replacement by SBM (SM25), and the fourth diet, fishmeal with 25% replacement by SBM and the inclusion of 2% of GroBiotic®-A (SM25-P). Each diet was supplied in three tanks (45 fish per diet); therefore, we had a 2 × 2 factorial design including two levels of FM replacement (0 and 25%) and two levels of prebiotics (0 and 2%), with three replicates (three tanks). The fishes were fed daily to apparent satiation, at each feeding period, all the uneaten feed in each tank was removed after one hour from feeding and was dry weighed (dw). Consumed feed was = total feed offered dw − uneaten feed dw, in that feeding period. The consumed feed was divided by the number of alive fish per diet.
The thermal growth coefficient (TGC) [50] and feed conversion rate (FCR = ingested food dry weight/fish wet weight gained) were estimated using initial and final weights for all fish and analyzed by two-way ANOVA with a statistical significance of P < 0.05, using STATISTICA v.8.0. After 56 days of the experiment, the fish were first anesthetized by immersion in seawater with tricaine methanesulphonate (MS-222) (90 mg/L) and then euthanized through brain puncture. The livers were extracted, embedded in RNA later, and incubated overnight at 4 °C, then they were stored at − 80 °C until used for RNA extraction.

RNA extraction and sequencing
For RNA sequencing, three biological replicates per diet were used for a total of twelve liver samples. Total RNA was extracted from 10 mg of every liver sample using the commercial kit Pure Link RNA mini kit (Invitrogen) following the manufacturer's instructions. The concentration and purity of the RNA were evaluated using the NanoDrop 2000 spectrophotometer (Thermo Scientific), the quantitation was verified using the Invitrogen Qubit fluorometer (Thermo Scientific). The integrity of the RNA was calculated with the 2100 Bioanalyzer Instrument (Agilent Technologies). The synthesis of complementary DNA (cDNA) to build the sequencing libraries was performed for each sample using the Illumina TruSeq RNA Sample Preparation Kit v2, according to the manufacturer's instructions. Paired-end RNA sequencing was performed using the Illumina Hiseq 4000 system to obtain reads of 2 × 100 nucleotides at the

Transcriptome assembly
The quality of the sequencing reads was analyzed using the FastQC software, then adaptors, ambiguous nucleotides, and low-quality sequences were eliminated with the Trimmomatic-0.36 software using a sliding window of 5 nucleotides and a minimum average quality of 28 in the Phred scale [51]. Finally, high-quality reads from all the libraries were assembled de novo using the Trinity v2.4.0 software with default parameters [52].

Differential gene expression analysis
Transcript abundance for each library was calculated following the two steps script (align_and_estimate_abundance) included in Trinity software: first Bowtie2 [53] was used for the alignment of the reads on the assembled transcriptome; then transcript quantification was carried out using the RSEM software [54]. A matrix of isoform counts was generated and used to analyze the differential expression among the libraries with DESeq2 [55] in R software [56] with filtering thresholds set to FDR below 0.01 and fold change above 2, following the run_DE_analysis script of Trinity software. The results of the differential expression analysis were represented with heatmaps which were generated following the analyze_diff_expr script included in Trinity software [52]. Finally, Venn diagrams were used to illustrate the differential expression results due to the effect of the different factors and their interaction using the software Venny 2.1.0 (http:// bioin fogp. cnb. csic. es/ tools/ venny/ index. html) [57].

Functional annotation
The functional annotation of the assembled transcriptome was carried out through the search of homologs with BLASTx [58] using the public databases UniProt (release 2019_2). The E value filter for this analysis was set at E-5.
Besides, the annotation with terms of gene ontology (GO) was made with the Blast2GO program [59]; each GO was classified into the categories of biological processes, cellular components, and molecular functions. The enrichment analysis for differentially expressed genes (DEGs) due to the effect of the factors (the FM replacement, prebiotics inclusion, and their interaction) was performed through Fisher's exact test (FDR < 0.05) using the Blast2GO program. As well, the metabolic pathways from the Kyoto Encyclopedia of Genes and Genomes database (KEGG) were downloaded through DAVID 6.8 [60], using DEGs. Additional enrichment analysis of biological processes and KEGG pathways was carried out for the whole liver transcriptome in DAVID 6.8 (P value < 0.01).

Growth performance
The weight gained, and the percent weight gained, resulted significantly higher due to the inclusion of prebiotics. Besides, the TGC was significantly higher in those fish fed with prebiotics. Additionally, no interaction was found between both factors (meal/prebiotic) and no significant differences were observed in the food consumed between treatments ( Table 2).

Sequencing and transcriptome assembly
A total of 712,373,692 high-quality paired reads of 100 bp were obtained including all libraries, and a total of 155,247 contigs were assembled with an average length of 972 bp, an N50 of 1722 bp, and a GC content of 47%, in total 150,952,398 bases were assembled.

Differential gene expression
The FM replacement caused the differential expression of 225 transcripts, the prebiotics inclusion caused differential expression on 242 transcripts, and the interaction of factors caused differential expression on 62 transcripts, for a total of 445 DEGs due to the factors and their interaction (Fig. 1). From these three categories, the top 50 DEGs were obtained, clustered, and represented in heatmaps. In the first heatmap, with DEGs corresponding to the effect of FM replacement, two main clusters were obtained: one grouping both diets with FM (FM and FM-P) and another grouping the diets including SBM (SM25 and SM25-P); replicate samples grouped congruently within every diet (Fig. 2). In the second heatmap, including DEGs corresponding to the effect of prebiotics inclusion, three main clusters were obtained: one grouping the control diet FM with SM25-P, the FM-P diet in an adjacent branch, and the SM25 diet as the outgroup. Therefore, the inclusion of prebiotics caused a similar differential expression pattern in the control (FM) and SM25-P diets. Again, all the replicate samples clustered congruently within the diets (Fig. 3). In the heatmap for the interaction of factors the FM and SM25-P diets were separated into different clusters and their replicate samples were clustered congruently; however, in this heatmap the replicates of the diets FM-P and SM25 did not cluster congruently with any factor (Fig. 4). The list of DEGs (up-and downregulated by each factor and their interaction) with annotations and relative expression (fold change) is provided as supplementary material (Supplementary Table 1).

Functional annotation
In the functional annotation through BLASTx, homologs were searched for the 155,247 contigs, and 59,027 matches were found with significant homology in the UniProt protein database, of which 51,609 had a GO and functional annotation. The species with the highest number of homologies was Homo sapiens (28,787), followed by Mus musculus (15,788) and Danio rerio (9183). Other fish species with significant matches were Takifugu rubripes (432), Oncorhynchus mykiss (388), and Salmo salar (279). Based on the results of functional annotation, GO terms of cellular component, molecular function, and biological processes were assigned. The most abundant terms in the cellular component category corresponded to intracellular part (38,744), intracellular organelle (32,680), and membrane-bounded organelle (29,361). In the category of molecular function, the top terms were protein binding (23,390), organic cyclic compound binding (15,094), and ion binding (14,913). Finally, in the category of biological processes, the most abundant terms were cellular metabolic process (28,175), metabolic process of nitrogen compounds (25,837), cellular component organization (15,651), and biosynthetic process (15,013). Important metabolic processes were enriched for the liver transcriptome like phospholipid biosynthetic process, steroid metabolic process, and protein transport, as well as the immune system process. KEGG pathways were also enriched including phosphatidylinositol signaling system, inositol phosphate metabolism, aldosterone synthesis and secretion, mTOR signaling pathway, B cell receptor signaling pathway, thyroid hormone signaling pathway, glucagon signaling pathway, and insulin signaling pathway among others (Table 3).

Enrichment analysis for DEGs affected by each factor
Using the DEGs obtained due to the effect of prebiotics inclusion, the biological processes with the best enrichment were Golgi vesicle transport, localization, tissue remodeling, and negative regulation of cell death. The DEGs obtained by the effect of the FM replacement enriched biological processes like ion transmembrane transport, response to stimulus, aromatic compound catabolic process, and actin filament organization. DEGs that were obtained by the interaction of factors enriched biological processes like immunoglobulin mediated immune response, protein autoubiquitination, single-organism catabolic process, and glucose metabolic process. The gene products included in every enriched process are shown in Table 4.

Discussion
In this study, we analyzed the liver gene expression of the Pacific yellowtail S. lalandi juveniles fed with diets containing different protein sources, with and without a prebiotic. The study was focused on the liver because it is the main Fig. 1 Venn diagrams representing the total number of DEGs (including isoforms) corresponding to the effects of A FM replacement (25%) by soybean meal, B prebiotics inclusion, and C the combination of factors in diets supplied to S. lalandi juveniles organ involved in the use of nutrients and is the center of intermediary metabolism in animals [33].
Regarding the experimental diets, the level of FM replacement by SBM tested in this study (25%) was high, considering that up to 20% is the recommended FM replacement level using SBM or soy protein concentrate for juveniles [61,62] and that inclusion of 30% or higher of SBM in diets have resulted detrimental for this species [63,64]. Although the level of FM replacement was high, it did not affect fish growth. Interestingly, the inclusion of prebiotics in diets with and without FM replacement improved the growth. Therefore, this replacement level using SBM when combined with GroBiotic®-A represents a good balance between cost and benefit for S. lalandi juveniles.
In the transcriptomic analysis, we distinguished DEGs potentially involved in metabolic pathways, growth, and the immune system. Prebiotics inclusion caused the major number of DEGs (242), followed by the FM replacement (225). On the other hand, the combination of factors (SM25-P diet) caused the least amount of DEGs. This suggests that prebiotics in SM25-P compensated the effect of SBM on the differential expression pattern, by keeping the expression of most of the genes at similar levels as those of the control diet.

Growth
Regarding the effect on growth, the diets including prebiotics (FM-P and SM25-P) resulted in the highest percentage of weight gain, the greatest weight gain in grams, and TGC. In these diets, we observed an upregulation on genes related to growth in fishes like the TGFBR3 [65,66] and EPS8 [67]. Moreover, the TGFBR3 gene has shown a significant association with the FCR in chickens [68]. On the other Fig. 2 Heatmap and clustering of DEGs corresponding to the effect of FM replacement (25%) by soybean meal in diets supplied to S. lalandi juveniles. Values in log2(TPM). Row labels indicate the UniProt IDs of the gene products and column labels indicate the name of the diets hand, growth was not affected by the FM replacement. In the SM25 diet, besides EPS8, the mTOR gene was upregulated, which is also related to fish growth [69]. The high expression of these genes and the activation of related metabolic pathways in diets containing SBM may explain why growth was not affected by this FM replacement.

Carbohydrate metabolism
The capability of the liver to maintain the balance between the intake and release of glucose plays an essential role in the maintenance of the homeostasis of glucose in the blood [70]. In this study, the glucose metabolism was mainly affected by the combination of FM replacement by SBM and prebiotics inclusion. Both diets with SBM showed high expression of the H6PD gene, and in the SM25 diet, the PFKFB4 gene was highly expressed. The enzyme 6-phosphofructo-2-kinase/fructose-2,6-bisphosphatase 4, which is encoded by PFKFB4 can reduce hepatic glucose production by increasing glycolysis and inhibiting gluconeogenesis, thereby lowering blood glucose [70]. Interestingly, this gene was associated with the adaptation of carnivorous fish to diets including plant meals [71]. Moreover, the inclusion of prebiotics increased the expression of G6PD and H6PD genes, both participating in glucose oxidation [72,73]. Besides, the protein phosphatase 1 regulatory subunit 3C (ppp1r3cb), which limits the breakdown of glycogen [74] was downregulated by the effect of prebiotics. Results suggest that the effect of the partial FM replacement combined with the inclusion of prebiotics enhanced glycolysis, and the energy supply obtained through this pathway was enough to reduce gluconeogenesis, improving protein weight gain.

Lipid metabolism
A better understanding of the physiological mechanisms of lipid metabolism of fishes is critical to replace fishmeal and fish oil in formulated diets, and if combined with gene expression studies would allow the identification of genes and potential markers for fish selection. According to our annotation and enrichment analysis, many genes involved in lipid metabolism were identified, like the ATP8A2 gene, which is upregulated in human adipocytes [75]. This gene participates in the asymmetric distribution and translocation of phospholipids and seems to be involved in the formation of vesicles and the uptake of lipid signaling molecules [76]. In our study, ATP8A2 was upregulated in diets including SBM (SM25-P and SM25), which can tell us about an effort by fish to maintain balance in the distribution of lipids due to SBM inclusion in the diet. High expression of ATP8A2 can be suggested as a molecular marker of fishes adapting to vegetable diets. Other important genes for lipid metabolism like LIPG and LPIN2 showed differential expression due to the interaction of factors. The LIPG gene, which participates in the hydrolysis of high-density lipoproteins [77] was downregulated in the SM25-P diet; by contrast, LPIN2, which plays an important role in the control of fatty acids metabolism [78], was upregulated. However, the highest expression of the LPIN2 was detected in the FM-P diet. In the intestinal mucosa of the sea bass (Dicentrarchus labrax), LPIN2 was upregulated in fish fed with diets including the highest percentage of FM [79]. In this study, this was observed only when prebiotics was added to the FM (FM-P).
Considering the mentioned above, both factors separately and their interaction induced changes in the expression levels of genes related to lipid metabolism in the liver. However, it is necessary to evaluate the expression and activity

Immune response
Finally, another relevant process affected mainly by the interaction of factors was the immune response. The immune system in fish is an important mechanism and the first line for protection against pathogens [80]. In this work, we found genes with potential participation in the immune response within the top DEGs. For example, the complement component C9 (C9) was found upregulated in SM25-P. This gene plays a key role in the innate and adaptive immune response by forming pores in the plasma membrane of target cells [81]. The expression of C9 was detected in the spleen, stomach, intestine, and head kidney, with the highest levels detected in the liver of the southern catfish (Silurus  meridionalis) [82]. Likewise, evaluations in Larimichthys crocea under a Vibrio alginolyticus challenge showed that two components of the terminal complement system: C7 and C9, were expressed in many tissues of adult healthy fish, with the highest levels detected in the liver [83]. Moreover, in SM25 and SM25-P the upregulation of the C3 gene was detected, which is involved in inflammatory processes and response to bacteria [84]. These results suggest that the SM25-P diet enhanced the immune response mediated by the gene products of C9 and C3.

Recommendations
Based on the overall results, both diets including prebiotics improved the growth rate and showed high expression of important metabolic genes in the liver. However, the FM-P diet may increase production costs. On the other hand, the 25% of FM replacement by SBM did not affect growth, and the inclusion of prebiotics in diets with SBM seems to improve the immune response. Moreover, the SM25-P diet compensated for the effect of the FM replacement on the liver at the transcriptomic level since the differential expression pattern in fish fed this diet resembled that of fish fed the control diet. Therefore, the replacement of 25% of FM by SBM combined with prebiotics may represent a good balance between cost and benefit for juveniles of this species. However, it is necessary to evaluate the response of fish for a longer period, long-term effects in growth, physiologic performance, immune response, and meal quality in adult organisms, and evaluate the effect of this diet on other organs of the digestive tract. Data availability All the RNA-Seq raw reads were deposited into the Sequencing Read Archive (SRA) of NCBI with the accession number SRR10211853 to SRR10211864. The BioProject ID of our data is PRJNA575250 and the BioSample accession is SAMN12816772.