Comparative transcriptomic analysis reveals an association of gibel carp fatty liver with ferroptosis pathway

Background Fatty liver has become a main problem that causes huge economic losses in many aquaculture modes. It is a common physiological or pathological phenomenon in aquaculture, but the causes and occurring mechanism are remaining enigmatic. Methods Each three liver samples from the control group of allogynogenetic gibel carp with normal liver and the overfeeding group with fatty liver were collected randomly for the detailed comparison of histological structure, lipid accumulation, transcriptomic profile, latent pathway identification analysis (LPIA), marker gene expression, and hepatocyte mitochondria analyses. Results Compared to normal liver, larger hepatocytes and more lipid accumulation were observed in fatty liver. Transcriptomic analysis between fatty liver and normal liver showed a totally different transcriptional trajectory. GO terms and KEGG pathways analyses revealed several enriched pathways in fatty liver, such as lipid biosynthesis, degradation accumulation, peroxidation, or metabolism and redox balance activities. LPIA identified an activated ferroptosis pathway in the fatty liver. qPCR analysis confirmed that gpx4, a negative regulator of ferroptosis, was significantly downregulated while the other three positively regulated marker genes, such as acsl4, tfr1 and gcl, were upregulated in fatty liver. Moreover, the hepatocytes of fatty liver had more condensed mitochondria and some of their outer membranes were almost ruptured. Conclusions We reveal an association between ferroptosis and fish fatty liver for the first time, suggesting that ferroptosis might be activated in liver fatty. Therefore, the current study provides a clue for future studies on fish fatty liver problems. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-07621-2.


Background
Fish fatty liver is a common physiological or pathological phenomenon in aquaculture. The causes are complex and not well-known, and mainly include imbalance nutrition diet, environmental pollutants, physiological factors and genetic mutation [1]. Fatty liver diseases have been found in most main farmed fish, and caused many problems, such as low feed efficiency, immune response, flesh and nutritional quality effects [1][2][3][4][5]. The utilization of artificially formulated feeds can bring nutritional, physiological and ecological effects to fish [2,[6][7][8][9]. However, exceeding nutrition, improper artificial formulated diets and overfeeding have led to a growing concern of liver fatty problems, such as hepatocyte enlargement, lipid accumulation, steatosis, fibrosis and necrosis [10,11].
Fatty liver can be mainly classified as "nutritional fatty liver" or "oxidative fatty liver" in aquaculture. Nutritional fatty liver, principally caused by imbalanced nutrition supply, is common in farmed fish and generally not a pathological symptom. It can be alleviated through adjusting diet formulation or feeding. If no effective strategy was carried on, it could be turned to the "oxidative fatty liver" or directly induced to hepatic fibrosis and necrosis, which causes irretrievable damages [1]. Oxidative stress might lead to oxidative fatty liver, which would arouse severe liver damage [1,12]. Overall, excess dietary energy intake and severe peroxidation can cause fish metabolic imbalance and then result in fatty livers. Some researchers suggest that many pathways, such as target-of-rapamycin complex 1 (Torc1), AMP-activated protein kinase (AMPK), transcription factor EB (TFEB), peroxisome proliferator activated receptor (PPAR), P53, nuclear erythroid 2-related factor 2 (Nrf2), c-jun Nterminal kinase (JNK), toll-like receptors (TLRs), myeloid differentiation primaryresponse protein 88 (Myd88), and nuclear factor κB (NF-κB) signaling pathways, might be related to fatty liver caused by high-fat/carbohydrate or over-nutrition diet [11][12][13][14][15][16][17]. For example, the decreased AMPK pathway can suppress autophagy and then worsen lipotoxicity in tilapia fatty liver [11], and its activation can reduce the expression of genes related to lipogenesis in rainbow trout liver [16]. However, some results seem to be controversial. For instance, the activation of Nrf2, JNK and TLRs-Myd88-NF-κB signaling pathways could lead to inflammation and worsen tilapia liver injury [12], while cAMP-JNK/NF-κBcaspase signaling pathway could protect the fatty liver tissues from more serious damage though regulating the hemostasis phosphorylation of JNK protein in Japanese seabass [17]. JNK and NF-κB signaling pathways might play a dual role in the fish fatty liver. Therefore, the mechanism of hepatic lipid accumulation and fish fatty liver are remaining enigmatic. In addition, too many affected pathways were identified and few studies had explored which key pathway was associated with fish fatty live.
Gibel carp (Carassius gibelio) is one of the most important aquaculture species in China [18][19][20][21][22][23] and the production yields in China have increased to 2,755,632 tons in 2019 [24]. In lotus-fish farming ponds, we found some individuals of gibel carp had fatty liver. To find out the cause, we first analyzed the liver histological structures and lipid accumulation of normal liver and fatty liver. Then, we conducted comparative transcriptomic analysis between, and identified a pathway of ferroptosis that was significantly activated in fatty liver. Finally, we confirmed that ferroptosis was activated in fatty liver by qPCR analysis and mitochondria morphological observation. Our current study establishes an association between ferroptosis and fish liver fatty, which provides a clue for future studies on fish liver fatty problems.

Morphological changes in fatty liver
The morphology of fatty liver showed obviously different from that of normal liver (Fig. 1a). The whole liver was more hypertrophic and the hepatocytes of fatty liver were more enlarged (Fig. 1b), showing less hepatocytes on an area of 10 μm square (Fig. 1c) (p < 0.01). And then, we performed oil red O staining (ORO) to compare the lipid accumulation between normal liver and fatty liver. As we expected, fatty liver showed more than 2 times lipid accumulation than that in normal liver (Fig. 1d-e) (p < 0.01).

Transcriptomic differences between normal liver and fatty liver
The transcriptomes of six liver samples in two groups were obtained using BGISEQ-500 Iillumina sequencing plat and each sample produced an average of 10.16 Gb clean bases. An average of 83.42 and 69.46% clean reads were mapping to the gibel carp' genome and gene sets (Table 1).
Finally, a total of 37,077 unigenes were obtained, which included 32,679 known genes and 4398 novel genes. Principal component analysis (PCA) showed groupings between normal liver and fatty (Fig. S1a). The correlation coefficient (R 2 ) is > 0.93 in group and < 0.83 between groups (Fig. S1b), suggesting the best-fitting regression line for the technical replicates according to standards and best practices for RNA-Seq. The results demonstrate that though there are individual differences, the two groups have a totally different transcriptional trajectory (Fig. S1) [25,26].

Enriched GO terms in fatty liver
A total of 3480 differentially expressed genes (DEGs) were assigned with the correction of p-values using FDR ≤ 0.001 (Table S1). Compared to normal, 1997 genes in fatty liver were upregulated and 1483 genes were downregulated (Fig. 2a). All up or downregulated DEGs were separately annotated into 2568 GO Gene Ontology (GO) terms and 3014 GO terms, among them, as visualized in Venn's diagrams, 964 and 1410 GO terms were just associated with up or downregulated DEGs (Fig. 2b). Sixteen and 20 GO terms were significantly enriched with the correction of q-value ≤0.05 ( Fig. 2 c-d).
According to the LPIA method, a total of 316 KEGG pathways and 962 GO terms with biological process classification, which two shared at least 1 DEG, were selected to conduct the pathway-pathway interaction network (Table S4). Four significant pathways were identified with LPIA_v_1.pl [34] in Perl with 1000 iterations (Table 2) After the weight Aij was calculated at random walk, a correlational network was conducted between pathways that included 43 enriched pathways in KEGG enrichment analysis and four significant latent pathways (marked with red point) in LPIA method (Table S5). It was clearly showed that four significant latent pathways connected with each other in the network (weight scores > 0.1) (Fig. 3b). The pathways "Porphyrin and chlorophyll metabolism" (23 DEGs, ko00860) ( Fig. S2a) and "HIF-1 signaling pathway" (49 DEGs, ko04066) ( Fig.  S2b) had more connections with other pathways in KEGG enrichment analysis. Among them, excepting ZIP8/14 & Ferritin with a up & down expression, glutathione peroxidase 4 (gpx4), a negative regulated gene in "Ferroptosis" [35][36][37], was downregulated, while other DEGs in "Ferroptosis" were all upregulated (Fig. 4a). The results suggest that the expression levels of genes in pathway "Ferroptosis" might be significantly different between normal liver and fatty liver.

Ferroptosis is activated in fatty liver
To validate the different expressions of genes in "Ferroptosis" between normal liver and fatty liver, four marker genes, gpx4 [35][36][37][38][39][40], acyl-CoA synthetase long-chain family member 4 (acsl4) [41,42], transferrin receptor 1 (tfr1) [43], and Glutamate-cysteine ligase (gcl) [44], were selected for qPCR analysis. Consistent with the transcriptome result, gpx4 was downregulated, while the others were all upregulated (Fig. 4b). Ferroptosis can also be characterized with the morphology of mitochondria [32,45]. To observe the morphological differences of mitochondria between fatty liver and normal liver, we performed transmission electron microscopy analysis. There was no significant cellular dysfunction, but compared to normal liver, the mitochondria densities in fatty liver were more condensed, and some of outer mitochondrial membrane had been ruptured and seemed like single-membrane (Fig. 5).
Taken together, ferroptosis was more sensitive in the fatty liver.

Discussion
In this study, we found morphological changes in fatty liver, such as hepatocyte enlargement and lipid accumulation (Fig. 1). GO and KEGG enrichment analysis showed that activities related to lipid biosynthesis, degradation, accumulation, peroxidation or metabolism  were more active in fatty liver (Fig. 2e & Fig. 3a). Importantly, a pathway of ferroptosis was significant different between normal liver and fatty liver that might be associated with lipid-related activities (Figs. 4 & 5), suggesting an association between ferroptosis and fish fatty liver. Based on the current data, we suggest that a significant pathway of ferroptosis might be associated with fish fatty liver. Ferroptosis is actually an iron-catalyzed-lipid peroxidation disorder, and is related with lipid peroxidation. Previous research in fish suggested a possible link between ferroptosis and lipid. For example, the ironcatalyzed lipid peroxidation was characterized in zebrafish and cultured shrimp [46,47]. Enzymes, such as glutathione peroxidase, could reduce the lipid peroxides in tilapia and other fish [48,49]. Oxidants could damage mitochondrial membrane permeability and electron transport chain integrity in zebrafish and grass carp [50][51][52]. However, ferroptosis in fish is not clear. Now, we establish an association between ferroptosis and fatty liver through comparative transcriptomic analysis, marker gene expression and mitochondrion morphology observation in fish.
Ferroptosis is a new form of regulated cell death that depends on iron-and lipid-based reactive oxygen species (ROS), and has been implicated in both normal and pathological physiology which involves in various biological contexts of diverse species increasingly and widely [32,35,38,53]. It is totally different from other reported forms of cell death, such as apoptosis, autophagy, necrosis and pyroptosis [40,53,54], and associated with various liver problems [38,55]. In fish, it potentially The relative expression of gpx4, acsl4, tfr1 and gcl in qPCR analysis. β-actin is used as the normalizer. Each bar represents mean ± standard deviation (SD) (n = 3). Asterisks stand for the significant differences between normal liver and fatty liver. (*: p < 0.05 and **: p < 0.01). NL: normal liver, FL: fatty liver was responsible for the lethality of zebrafish VitEdeficient embryos [56]. It might be activated by low temperature (Nile tilapia), hypoxia (Silver sillago) stress and heavy metal (Japanese flounder) [57][58][59], and could be induced by Escherichia coli and then led to the red blood cells (RBCs) death of grass carp [60]. It could be checked at the gene expression level and characterized by the morphology of mitochondria [32,35,36,38,39,[41][42][43][44][45]61]. Some genes in ferroptosis are classified as markers. For example, gpx4 is required for the clearance of lipid ROS. If it was inactivated or inhibited, lipid ROS would accumulate, and then induced ferroptosis [35][36][37][38][39][40]. acsl4, as an essential component for ferroptosis execution, could dictate ferroptosis sensitivity. The more upregulated expression of acsl4, the more sensitive to ferroptosis [41,42]. In this study, compared to normal liver, the relative expression of negative regulator gpx4 was significantly downregulated and the relative expression of three positive actors (acsl4, tfr1 and gcl) were significantly upregulated, implying that fatty liver was more sensitive to ferroptosis. In addition, ferroptosis showed smaller, condensed and outer membrane ruptured of mitochondria with crista diminished or vanished in cellular morphological characteristics [62,63]. Moreover, the GO and KEGG enrichment analysis suggest lipidrelated activities were active in fatty liver. Therefore, our results indicate that ferroptosis might be associate with liver fatty and it is different between normal liver and fatty liver.
The main method used in this study was LPIA, which was more appropriate to our study in one challenge with metabolic relevant. It could give a direct understanding to the key cellular mechanisms in biological activities [34,64]. In this study, four latent pathways were identified by LPIA (Table 2), including ferroptosis. Among them, the remaining three latent pathways might also connect with ferroptosis (Fig. 3b). For example, "Porphyrin and chlorophyll metabolism" pathway was associated with the energy, ion and erythroid heme synthesis regulation, especially in responding to various stresses [65, . "HIF-1 signaling" pathway was involved in cellular responses to low oxygen environments and growth factors [67,68]. It was noteworthy that both "Porphyrin and chlorophyll metabolism" and "HIF-1 signaling" pathways had more relations with other enriched pathways (Fig. 3b). Minerals were fundamental nutrients and mineral absorption may depend on dietary ingredient composition. Some genes in "mineral absorption pathway" were involved in iron metabolism [69,70]. These could explain why they were all identified and suggest that ferroptosis may be the key pathway involved in fatty liver.
Overall, as lipid accumulation continued to increase, lipid peroxidation would also increase [52]. If hepatocytes could not be sufficient capacity to eliminate lipid peroxides, ferroptosis might be activated [32,53]. Certainly, more research on the mechanism underlying ferroptosis and fish fatty liver need to be carried out. Since ferroptosis could be regulated and numbers of small molecule inhibitors, such as Vitamin E, the nature's most-efficient ferroptosis inhibitor, have been identified [32,33,40,56,[71][72][73], which suggests new prevention strategies and promising therapies of fish fatty liver.

Conclusion
Based on detailed comparison of histological structure, lipid accumulation, transcriptomic profile, ferroptosis marker gene expression, and hepatocyte mitochondria between normal liver and fatty liver, our study reveals an association between ferroptosis and fatty liver in fish for the first time, suggesting that ferroptosis might be activated in fish liver fatty. The current study provides a clue for future studies on fish fatty liver and effective prevention strategies in fish fatty liver problems.

Animals and sample collection
All experimental procedures in this study were performed in accordance with the guidelines and after approval of the Animal Care and Use Committee of Institute of Hydrobiology, Chinese Academy of Sciences (IHB, CAS, Protocol No. 2016-018). A total of 3000 juveniles of allogynogenetic gibel carp, which came from same clones, were transported from National Aquatic Biological Resource Center of Institute of IHB, CAS, which is located in Wuhan, China, to the Luo-Fu Lotus Farm, which is located in Jinggangshan, China, for lotusfish culturing. Fifty juveniles were selected randomly to weight (Table S6). Half of the juveniles was set as an overfeeding group, which was continuing feeding until without fish activities in apparent satiation two times per day. The other half was set as a control group, which was cultured without artificial feeding. All the juveniles were cultured in lotus-fish culturing ponds at a density of 150 juveniles per 667 m 2 . Water samples were collected from three randomly ponds three times (April 22th, June 3th, October 12th, 2019). Eight physicochemical water quality parameters (water temperature, dissolved oxygen, pH, total nitrogen, ammonia, nitrate, nitrite and total phosphorus) were analyzed (Table S7). After 180 days, commercial fish were harvested. Fifty commercial fish of each group were selected randomly to weight (Table S6). Compared to the individuals in the control group with normal liver, most of individuals in the overfeeding group had bulge belly and randomly selected of them had fatty liver in general appearance after dissecting. Then, the liver samples of three individuals of the control group with normal liver and three individuals of the overfeeding group with fatty liver were collected randomly for histopathological observation and molecular analysis.
All fish were first totally immersed and deeply anesthetized with tricaine methanesulfonate (MS-222, 35-40 mg/L, Servivebio, Wuhan) until losing all rhythmic opercular movements for a minimum of 30 min, then sacrificed by rapidly cutting off the spinal cord adjacent to the head. Procedures were performed to minimize fish suffering as far as possible. All sections of this study were carried out in compliance with the ARRIVE guidelines [74], and a completed ARRIVE guidelines checklist was included (Checklist S1).

Histological structure, lipid accumulation and morphological observation
Liver tissues were cut in sections and fixed in 4% paraformaldehyde (PFA) overnight at 4°C, and subsequently embedded in paraffin. One part of sections was stained with hematoxylin & eosin (HE) (Beyotime, Suzhou) and performed as previously described [75,76]. The Carl-Zeiss microscopy (Analytical & Testing Center, IHB, CAS) was used for histological structure observation and photomicrographs.
The other sections were used for detecting the lipid droplet morphology by ORO staining according to the previous method [77,78]. The sections were totally dehydrated with 30% sucrose, and then embedded in optimum cutting temperature (OCT) compound. After being rapidly frozen and cut in Thermo CRYOSTAR NX50 microtome at 10 μm thickness, the sections were dyed with ORO (Thermo, USA) and counterstained with hematoxylin (Servivebio, Wuhan). Via the amount of ORO staining from 3 frames per biopsy, liver lipid accumulation was quantified using pixel numbers with Image J [79].
The specimens for morphological observation by electron microscopy were prepared as described previously [80,81]. Liver tissues were fixed with 2.5% glutaraldehyde for 24 h at 4°C, and in 1% osmium tetroxide (OsO4) for 2 h at 4°C, and then gradiently dehydrated with ethanol. After that, sections were embedded in epoxy resin Epon812 for overnight and cut in Leica DMIRB ultrathin microtome at 60-80 nm thickness, stained with uranyl acetate and lead citrate, and observed with a HC-1 80.0 KV Hitachi TEM system (Analytical & Testing Center, IHB, CAS).
RNA extraction, sequencing, assembly and functional annotation of liver transcriptomes Three liver samples of two groups, normal liver and fatty liver, were listed as biological replicates (NL-1, NL-2, NL-3 and FL-1, FL-2, FL-3). Total RNA was extracted with RNeasy Mini Kit (Qiagen, Beijing) according to the manufacture's protocols. Proper quality and quantity were checked with Nanodrop® ND-2000 spectrophotometer (LabTech, USA) and Technologies 2100 Bioanalyzer (Agilent Tech, USA) by measuring the 260/280 nm absorbance ratio. Fifty μg total RNA of each samples were used for cDNA library establishing and high-throughput sequencing via BGISEQ-500 platform in BGI Genomics Co., Ltd., Shenzhen, China. Clean reads were mapped to the gibel carp' genome and genes using HISAT (v2.1.0) [82] and Bowtie2(v2.2.5) [83] (unpublished data). The expression levels of transcripts were calculated with RSEM [84]. Pearson's correlation coefficient (r) and PCA were measured by using cor & princomp in R for the strength of the association between the two groups. DEG analysis was conducted by DEseq2 [85,86] with q-value (adjusted pvalue) ≤ 0.05. GO and KEGG pathway enrichment analysis were conducted with phyper in R. p-value was corrected with false discovery rate (FDR) cut-off of 0.01 and q-value ≤0.05 was used as the threshold to judge the significance of GO or KEGG enrichment. The raw sequences were deposited into the NCBI Sequence Read Archive (SRA) database (Accessions PRJNA675741).

Latent pathways network construction and the significant pathways identification
Latent pathway identification analysis, LPIA, developed by Pham et al., is a method to identify the interactions of pathways associated with DEGs and biological processes [34,64]. After preparing the three distinct but interrelated sources of biological pathways (such as "KEGG", P), biological functions (such as "GO", G) and gene transcriptional response in different conditions (such as "DEGs", DE) (Table S4), 1) First, we constructed a P-G bipartite graph if P and G shares a non-empty intersection of DE, we calculated the weighted edge, denoted W GP , which stands for the intersection of P and G as follow 2) And then, converted the P-G bipartite graph to P-P pathway network if two pathways, P i and P j shares non-empty intersection of G, the weight between P i and P j, denoted A i j is 3) Finally, measured pathway importance via eigenvector centrality to determine the significance of a node in the above network after proper iterative operation. The weight scores > 0.1 were used to draw P-P pathway network with cytoscape 3.7.2 [87]. The whole process was repeated using the bootstrap method [88] and the adjusted p-values were calculated to account for multiple testing with methods described by Dudoit and van der Laan [89].

Quantitative reverse transcription PCR (qPCR)
The first-strand cDNA was synthesized from total RNA following the protocol of Thermo Scientific™ RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, USA) in a 20 μl reaction volume. The expression level was analyzed on a CFX96™ Real-Time PCR System (Bio-Rad, USA) using an iTaqTM Universal SYBR®Green Supermix (Bio-Rad, USA). Gene-specific primers (Table  S8) were designed with Primer premier 5.0 [90]. The reaction system, protocol, and endogenous control selection were conducted as previously described [91][92][93]. All liver samples were analyzed with three biological replicates and the relative expression levels of target genes were normalized to β-actin and calculated by the 2 -ΔΔCT method. p-value < 0.05 was considered statistically significant. Additional file 2: Table S1 DEGs (FDR ≤ 0.001) in normal liver and fatty liver. Lists of DEGs include Gene ID, Length, FPKM, Log 2 fold change, Qvalue, Pvalue, and Annotation.
Additional file 3: Table S2 GO terms with all up/downregulated DEGs enrichment analysis. GO terms with all upregulated DEGs (GO terms with all up-DEGs), GO enrichment analysis with all upregulated DEGs (enriched GO-up), GO terms with all downregulated DEGs (GO terms with all down-DEGs) and GO enrichment analysis with all downregulated DEGs (enriched GO-down). NL: normal liver, FL: fatty liver.
Additional file 4: Table S3 KEGG pathway enrichment analysis. Pathway ID, Pathway Name, KEGG function classification, number of Candidate and total DEGs/genes, Rich Ratio, P-value, Q-value are shown. NL: normal liver, FL: fatty liver.
Additional file 11: Table S8 Primers used in this study.

Acknowledgments
We would like to thank Dr. Xiao Yuan at The Analysis and Testing Center of Institute of Hydrobiology, Chinese Academy of Sciences for her assistance with electron microscope analysis. The research was supported by the Wuhan Branch, Supercomputing Center, Chinese Academy of Sciences, China.
Authors' contributions XJZ performed laboratory experiments, analyzed the data and drafted the manuscript. LZ designed the studies and revised manuscript. WJL, WXD, XYL, ZWW and ZL participated in the laboratory experiments. YW participated in the data analysis. XYM and MD participated in the aquaculture conditions analysis. JFG conceived the studies and revised the manuscript. All authors have read and approved the manuscript.

Availability of data and materials
The raw sequences supporting the conclusions of this article were deposited into the NCBI Sequence Read Archive (SRA) database under the accession