The fundamental feed formula is shown in Table 1. Fish meal was used as main protein sources. VA acetate (500,000 IU/g, Maclin Biochemical Co., Ltd., Shanghai, China) is the supplemental source of VA. Three isonitrogenous and isolipidic diets were formulated adding VA 0, 200 and 16000 IU/kg (3 201, 4 730 and 16890 IU/kg, actually) denoted as LA, MA and HA groups respectively. All raw materials were screened for 80 mesh and thoroughly mixed. After adding 16% moisture, the pellets were made by the pelletizer at 120°C. After air drying, the feed were stored at −20°C.

Concentrated sulfuric acid, boric acid, hydrochloric acid, sodium hydroxide, petroleum ether, methanol, acetonitrile acetone, sodium sulfate and anhydrous ethanol purchased from Nanchang Jinsha chemical raw materials Co. LTD. Trizol reagent, electrophoretic buffer, Reverse transcription kit, loading buffer and SYBR Premix Ex Taq kite purchased from Jiangxi Guanyin Biotechnology Co. LTD.

The juvenile Monopterus albus were artificially bred by our team. All fish were reared temporarily for 14 days to acclimate to the conditions. After fasting for1 day, 480 fish (7.3 ± 0.02 g) were randomly assigned to 12 separate tanks of 3 groups (HA, MA, and LA, four repetitions per group) with 40 fish per tank. Each tank was filled with the right amount of fresh grass, Eichhornia crassipes. Fish were fed 3–5% of the body weight once daily at 18:00 p.m. A tuck net removed the feces and residual feeds. The breeding experiment lasted for 10 weeks and replaced with clean water weekly. The environmental conditions were: water depth 0.4 m; temperature, 28°C ± 1°C; dissolved oxygen >5.0 mg/L.

After the feeding trial, the fish were anaesthetized using MS222. Three fish from each tank were sampled randomly. The liver from twelve fish each group were dissected into 1.5 mL tubes and stored in −80°C for further RNA extraction. Three biological replicates were used for transcriptome sequencing and six biological replicates were used for metabolome detection. The remaining three fish each group were used to detect body composition.

The crude protein, crude lipid, moisture and ash in whole fish body were measured following the methods of the Association of Official Analytical Chemists. The crude protein was determined by Kjeldahl method (FOSS, Kjeltec TM 8200). The crude lipid was estimated by using Soxhlet extraction method (FOSS, Soxtec 2050). The moisture was determined by drying samples in an oven at 105°C until constant weight. The ash was determined by combustion in a muffle furnace at 550°C for 4 h.

Liver RNA was extracted using Trizol reagent (Invitrogen, CA, United States) according to manufacturer’s instructions. Agarose gel electrophoresis was used to monitor RNA degradation. The Nano-Drop 2000 spectrophotometer was used to detect RNA concentration and purity. Further, the Agilent 2100 RNA 6000 Nano kitwas used to assess RNA integrity. Higher quality RNA was selected for high-throughput sequencing. Three RNA replicates of each sample were used to build RNA-seq library. After total RNA was extracted, the first strand cDNA was synthesized by random hexers and reverse transcriptase, and the second strand cDNA was synthesized by DNA polymerase I and RNase H. Then, a cDNA library was constructed for sequencing. The cDNA library was sequenced using Illumina HiSeqTM4000 by Gene Denovo Biotechnology Co. (Guangzhou, China). The transcriptome results were compared pair-to-pair.

Six genes (Table 2) were selected randomly for validation of RNA-Seq data by qPCR using a SYBR Premix Ex Taq kit (Invitrogen) according to the manufacturer’s instructions. The same RNA samples were used for both Illumina library synthesis and the qPCR verification assay. The first strand cDNA was obtained from 2 μg of total RNA using a PrimeScript first strand cDNA synthesis kit (Takara, Dalian, China). Melting curve analyses were performed following amplification. The specific primers used for qPCR are listed in Additional Table 2, and 18 s gene was used as an endogenous control. The thermal profile for SYBR Green qPCR was 95°C for 90 s, followed by 40 cycles of 95°C for 5 s, 60°C for 15 s, 72°C for 20 s.

The Liver samples were thawed at 4°C, appropriate samples were added into the pre-cooled methanol/acetonitrile/water solution (2: 2, v/v), vortex mixing, ultrasonic for 30 min, standing at −20°C for 10 min. After 14,000 g at 4°C for 20 min, supernatant vacuum drying, adding 100 μL acetonitrile solution during mass spectrometry (acetonitrile: Water = 1:1, v/v) redissolved, vortex, 14,000 g centrifuged at 4°C for 15 min. Then, supernatant was taken for analysis. Analysis was performed using an UHPLC (1290 Infinity LC, Agilent Technologies) coupled to a quadrupole time-of-flight (AB Sciex TripleTOF 6600) in Shanghai Applied Protein Technology Co.,Ltd. The metabolome results were compared in pairs.

Firstly, differential genes and differential metabolites related to lipid metabolism were identified by KEGG pathway analysis. Secondly, O2PLS (two-way orthogonal projections to late structures) model was established to analyze the relationship between them. Finally, pearson correlation coefficients were used to assess the correlation between differential genes and differential metabolites.

The results were expressed as mean ± SD. (standard deviation of the mean) and were analyzed using one-way ANOVA in SPSS 19.0 software. Any significant differences were further investigated by comparing the group means using Turkey’s test. Statistical significance was considered if p < 0.05.

The crude lipid content in the whole body was significantly different (p < 0.05) among the different treatments, and that in the whole body of the LA group was significantly lower than that of the MA and HA groups (Table 3). The crude lipid content of the whole body increased with increasing dietary VA levels. In contrast, the differences in crude protein, moisture, and ash in the whole body were not significant.

To validate the transcriptome data, 6 DEGs were selected randomly for qPCR. Althought the test DEGs displayed different expression levels (Figure 1), in general, the qPCR results showed a positive correlation with transcriptome, indicating the reliability and accuracy of the transcriptome analysis.

After removing the low-quality reads, 50.52 million clean reads were obtained. A total of 32.0 Gb of data were generated by Illumina sequencing. The results of transcriptome have been submitted to NCBI (SRA accession no. SRP398570).Comparative transcriptome analysis identified 1323 differentially expressed genes (DEGs) among the three groups (Figure 2). Based on the KEGG annotation, lipid metabolism-related DEGs and pathways were further analyzed. A total of 46 DEGs were enriched in 14 lipid metabolism-related pathways (Table 4). Among them, LAvsMA contained 10 pathways and 16 DEGs; MAvsHA contained 14 pathways and 36 DEGs; LAvsHA contained 14 pathways and 46 DEGs. Significance analysis revealed that sphingolipid metabolism (ko00600), glycerol metabolism (ko00561), steroid hormone biosynthesis (ko00140), and primary bile acid biosynthesis (ko00120) were significantly enriched. There were 16 DEGs in the four significantly enriched pathways. Due to the crude lipid content was continuously rising, further trend analysis of the screened DEGs was necessary. Trend analysis revealed that PLPP1a, PLPP2b, and CYP21a2 were continuously upregulated, AKR1D1, UGT1a1, and CYP1a were continuously downregulated. The relationships between six genes and the four significantly enriched pathways are shown in Table 5.

According to analyses of the metabolites, 10878 metabolites and 8540 metabolites were obtained in positive and negative modes, respectively. Based on the variable important in projection (VIP) values of OPLS-DA, the significantly different metabolites were screened based on a VIP ≥1 and p ≤ 0.05. In negative ion modes, a total of 201 different metabolites were detected in the three groups (Figure 3). To explore the lipid metabolism pathways, we selected different lipid-relevant metabolites for KEGG enrichment analysis. The five lipid metabolism pathways were mainly related to arachidonic acid metabolism, primary bile acid biosynthesis, unsaturated fatty acid biosynthesis, sphingolipid metabolism, and steroid hormone biosynthesis. There were six metabolites identified in the five pathways: eicosapentaenoic acid, sphinganine, 11-beta-hydroxyprogesterone, hydroxyeicosatetraenoic acid, cholic acid, and glycochenodeoxycholate (Table 6). In positive ion modes, no metabolites were detected.

To screen for DEGs and different metabolites associated with lipid metabolism, integrated transcriptome and metabolome analyses were performed using the O2PLS, correlation coefficient models and KEGG pathway. The common KEGG pathways were analyzed for the differential genes and metabolites. Integrative analysis revealed that primary bile acid and steroid hormone biosynthesis were enriched in both the transcriptome and metabolome.

Transcriptomic and metabolomic data were correlated analysis by O2PLS model. Components of the O2PLS model were calculated and reflected the contribution of the transcriptome and metabolome to the total variation. Therefore, the O2PLS model was very reliable for explaining the total variation in the comprehensive analysis of transcriptome and metabolome.The combined loading diagram of metabolites and transcripts was constructed according to the loading values shown in Figure 4. High absolute loading values indicated a strong correlation between the differential genes and differential metabolites. DEGs and differential metabolites farther from the origin were more closely related to each other. The correlation between metabolites and genes were evaluated by calculating pearson coefficients. The correlation matrix of the heat map revealed negative (blue) and positive (red) correlations between lipid-related differential metabolites and genes (Figure 5). Correlation network maps showed lipid-related differential genes or metabolites at important associated positions (Figure 6). Sphinganine was the main associated metabolite and it was correlated with five DEGs.

The primary bile acid biosynthesis pathway was significantly enriched in both the transcriptome and metabolome, indicating that there is a complex interaction between VA and bile acids. Retinoic acid receptor (RAR), farnesoid X receptor (FXR), and retinoid X receptor (RXR) are the main regulators of VA and bile acid homeostasis (Zoubek et al., 2017). In the liver, VA and its metabolites activate RAR and RXR, generate RAR/RXR heterodimers, and regulate bile acid synthesis and metabolism (Li et al., 2021). In addition, retinoic acid can inhibit the genes involved 5in bile acid synthesis (Mamoon et al., 2014), in which RAR and RXR may be involved. Studies have shown that VA inhibits the expression of the bile salt synthesis rate-limiting enzyme CYP7A1, thereby reducing bile generation (Jahn et al., 2016). In this study, AKR1D1 was inhibited instead of CYP7A1. AKR1D1 plays an important role in steroid metabolism (Gathercole et al., 2021), but biologically, it is better known for its involvement in bile acid biosynthesis, where it catalyzes a key step that introduces the 5β-configuration into primary bile acids (Russell, 2009). This decrease in cholic acid levels may be due to the inhibition of AKR1D1. Another reason is that cholic acid is used to produce 12α-hydroxylated bile acid, which lead to hepatic lipid accumulation (Maegawa et al., 2021). In contrast, the glycochenodeoxycholate content was found to be increased in the metabolome. Glycochenodeoxycholate is a hydrophobic bile salt and that can cause lipid accumulation in the liver (Chen and Chen, 2014). Simultaneously, glycochenodeoxycholate can increase lipid solubility and promote intestinal absorption (Wu et al., 2016). Notably, high concentrations of glycochenodeoxycholate induced hepatocyte apoptosis (Wang et al., 2005). This may also explain why high levels of liver fat are disadvantageous in fish.

Fatty acids are components of various molecules, including phospholipids, sphingolipids, and esters. The regulation of fatty acid metabolism by VA has been confirmed in other studies (Chen and Chen, 2014; Blaner, 2019). The genes for fatty acid extension, desaturation and synthesis pathways can be regulated by VA. This is achieved through retinoic acid (RA), a product of VA metabolism.

The first step in fatty acid synthesis is the conversion of acetyl-CoA to malonyl-CoA by acetyl-CoA carboxylase. It is then catalyzed by fatty acid synthase to form palmitic acid. Palmitic acid can be prolonged and desaturated to produce fatty acids with longer chains and double bonds (Kihara, 2012; Zhang et al., 2016). Elongases add two carbon atoms each time to create longer-chain fatty acids. Desaturases introduce double bonds into saturated and unsaturated fatty acids, resulting in the production of monounsaturated and polyunsaturated fatty acids, respectively. The mechanism of VA regulate fatty acid synthesis has been extensively studied. Studies have shown that when VA deficiency, the mRNA levels of acetyl-CoA carboxylase, fatty acid synthase and stearoyl-CoA desaturase 1 were reduced; when RA treatment, the mRNA levels were increased (Yang et al., 2021).

Ayuso et al. (2015) found that VA increases the content of polyunsaturated fatty acids in pork. Similar to the results of our study, the levels of the polyunsaturated fatty acids, eicosapentaenoic and hydroxyeicosatetraenoic acid, increased significantly. This is probably because of the increase in stearoyl-CoA desaturase 1. Stearoyl-CoA desaturase 1 was responsible for the formation of double bond and regulation of lipogenesis in the process of fatty acid synthesis (Koeberle et al., 2016; ALJohani et al., 2017). RA could induce stearoyl-CoA desaturase 1 mRNA expression (Miller et al., 1997) and desaturase index in adipose tissues (Gorocica-Buenfil et al., 2008).

Oliveros et al. (2007) found that VA deficiency induced a hypolipidemic effect by decreasing the total phospholipid content in the liver, owing to low phosphatidylcholine synthesis and enhanced fatty acid oxidation in rats, and this was consistent with the results of our study. In the present study, PLPP1a and PLPP2b were significantly upregulated with increasing VA content. PLPP1a and PLPP2b are phospholipid phosphatases that belong to the lipid phosphate phosphatase family. Lipid phosphate phosphatase has wide substrate selectivity and can hydrolyze a variety of lipid phosphoric acids containing single lipid bonds, including phosphatidic acid, lysophosphatidic acid, diacylglycerol pyrophosphate, and sphingosine 1-phosphate, and it participates in the synthesis of glycerol (Toke et al., 1998; Li et al., 2017; Tang and Brindley, 2020). The significant increase in sphinganine levels in the metabolome may be related to the hydrolysis of sphingosine 1-phosphate.

VA is essential for the synthesis and secretion of steroid hormones. The expression of 11β-hydroxysteroid dehydrogenase type 1 was increased by VA deficiency in the hypothalamus (C62.5%) and hippocampus (C104.7%) (Marissal-Arvy et al., 2013). RA treatment could reverse the effects of VA deficiency on 11β-hydroxysteroid expression. The inhibitory effect of retinoids on 11β-hydroxysteroid dehydrogenase type 1 expression has been demonstrated peripherally in the liver and visceral fat of WNIN/Ob obese rats (Sakamuri et al., 2011), in the adipose tissue of mice (Sakuta et al., 2006), and in myotubes in vitro (Aubry and Odermatt, 2009). In contrast, in this study, the level of 11β-hydroxyprogesterone was increased by VA, but this discrepancy may be species related. Within the cell, 11β-hydroxysteroid dehydrogenase type 1 catalyzes the regeneration of active glucocorticoids, thereby amplifying their action (Marissal-Arvy et al., 2013). Similarly, in our study, we found that CYP21a2 was associated with glucocorticoid synthesis (Heidelberg, 2008) which was significantly upregulated.

In this study, comparative metabolomics and transcriptomics were applied to investigate lipid metabolite changes in rice field eel when fed different levels of VA. Comparative transcriptome analysis revealed that VA affected multiple lipid metabolic pathways, especially sphingolipid and glycerin metabolism, and steroid hormone and primary bile acid biosynthesis. Metabolite analysis revealed that the levels of lipids, such as eicosapentaenoic acid, sphinganine, 11-beta-hydroxyprogesterone, hydroxyeicosatetraenoic acid, and glycochenodeoxycholate, increased significantly. Integrative analyses of the transcriptome and metabolome revealed a significant correlation between differential metabolites and differential genes, suggesting that these are the main metabolic pathways. These findings provide reliable models at the transcriptional and metabolic levels which will facilitate further investigation into the effects of VA on lipid metabolism.