Integration of transcriptomics and metabolomics reveals the molecular mechanisms underlying the effect of nafamostat mesylate on rhabdomyolysis-induced acute kidney injury

Objective: To investigate the role and mechanisms of action of nafamostat mesylate (NM) in rhabdomyolysis-induced acute kidney injury (RIAKI). Methods: RIAKI rats were assigned into control group (CN), RIAKI group (RM), and NM intervention group (NM). Inflammatory cytokines and proenkephalin a 119–159 (PENKID) were assessed. Cell apoptosis and glutathione peroxidase-4 (GPX4) were detected using TUNEL assay and immunohistochemical staining. Mitochondrial membrane potential (MMP) was detected by JC-1 dye. The expression of genes and metabolites after NM intervention was profiled using transcriptomic and metabolomic analysis. The differentially expressed genes (DEGs) were validated using qPCR. The KEGG and conjoint analysis of transcriptome and metabolome were used to analyze the enriched pathways and differential metabolites. The transcription factors were identified based on the animal TFDB 3.0 database. Results: Serum creatinine, blood urea nitrogen, and PENKID were remarkably higher in the RM group and lower in the NM group compared to the CN group. Pro-inflammatory cytokines increased in the RM group and notably decreased following NM treatment compared to the CN group. Tubular pathological damages were markedly attenuated and renal cell apoptosis was reduced significantly in the NM group compared to the RM group. The expression of GPX4 was lower in the RM group compared to the CN group, and it increased significantly after NM treatment. A total of 294 DEGs were identified in the RM group compared with the NM group, of which 192 signaling pathways were enriched, and glutathione metabolism, IL-17 signaling, and ferroptosis-related pathways were the top-ranking pathways. The transcriptional levels of Anpep, Gclc, Ggt1, Mgst2, Cxcl13, Rgn, and Akr1c1 were significantly different between the NM and RM group. Gclc was the key gene contributing to NM-mediated renal protection in RIAKI. Five hundred and five DEGs were annotated. Compared with the RM group, most of the upregulated DEGs in the NM group belonged to Glutathione metabolism, whereas most of the downregulated DEGs were related to the transcription factor Cytokine-cytokine receptor interaction. Conclusion: NM protects the kidneys against RIAKI, which is mainly associated with NM mediated regulation of glutathione metabolism, inflammatory response, ferroptosis-related pathways, and the related key DEGs. Targeting these DEGs might emerge as a potential molecular therapy for RIAKI.


Introduction
Rhabdomyolysis (RM) is a clinical syndrome of a muscular disorder caused by trauma, stroke, major artery occlusion, status epilepticus, genetic defects, infections etc. Rhabdomyolysis is characterized by the release of skeletal muscle-cell contents, such as myoglobin, sarcoplasmic proteins, and electrolytes into circulation (Shapiro et al., 2012;McMahon et al., 2013). Elevated myoglobin indicates the severity of RM (Bagley et al., 2007), which can result in life-threatening complications, including acute kidney injury (AKI), hyperkalemia, hypocalcemia or hypercalcemia, metabolic acidosis, hyperuricemia, hyponatremia, hypovolemic shock, and cardiac dysrhythmias (Bosch et al., 2009;Cervellin et al., 2010). It is worth noting that RM induced AKI (RIAKI) occurs in 19%-58% of patients with rhabdomyolysis and is associated with high mortality (Simpson et al., 2016).
The pathophysiology of RIAKI is thought to be associated with the elevation of circulating myoglobin. Under normal physiological conditions, myoglobin binds to plasma globulins and only a very small portion of myoglobin is released into circulation and reaches the glomeruli. However, in severe muscle damage, the circulating myoglobin overwhelms the globulins' binding capacity, which enables the free myoglobin to reach the renal glomeruli and tubules (Krouzecky et al., 2003;Petejova et al., 2014). The presence of myoglobin in the urine (known as myoglobinuria), leads to tubular cast formation and intratubular obstruction, accumulation of iron, and proximal tubular cell injury. It has been shown that myoglobin can bind to the Tamm-Horsfall protein in the tubular cast, particularly in an acidic environment. The cast-bound form of myoglobin results in oxidative stress, which leads to vasoconstriction and a significant reduction of renal cortical blood flow (Moore et al., 1998;Giannoglou et al., 2007). Another factor of pathogenesis of RIAKI is volume depletion caused by sequestration of fluids distributed into injured muscles. Cellular constituents released from damaged muscles results in metabolic acidosis and the resulting aciduria further aggravates tubular injury due to low urinary pH and presence of globin and ferrihemate derived from myoglobin can cause a direct nephrotoxic effect (Chatzizisis et al., 2008).
Currently, there is no specific treatment for RIAKI. Dialysis remains the core treatment for patients with AKI secondary to rhabdomyolysis. It is hypothesized that dialysis can remove the circulating myoglobin released from injured muscles and metabolic toxic substances in the circulation. Continuous renal replacement therapy (CRRT) has emerged as one of the most important modalities for the prevention and treatment of RIAKI, especially for hemodynamically unstable patients (Vanholder et al., 2011). Due to the rapid repletion of plasma myoglobin following a single hemodialysis therapy, intermittent hemodialysis is usually ineffective in maintaining a low level of plasma myoglobin, while CRRT has greater clearance of myoglobin, and is more effective than conventional filters (Vanholder et al., 2011;Zimmerman and Shen, 2013). To ensure the safety and efficacy of CRRT, prolonged anticoagulation therapy is required to prevent clotting in the extracorporeal circuits. Heparin has been used frequently as an anticoagulant in the implementation of CRRT, but may increase the risk of hemorrhage and thrombocytopenia (Webb et al., 1995;van de Wetering et al., 1996). Nafamostat mesylate (NM) is a synthetic serine protease inhibitor, which exerts an inhibitory effect on various enzyme systems, such as the fibrinolytic, coagulation, and complement system (Maruyama et al., 1996). In 1990, NM was introduced into CRRT as an alternative anticoagulant in patients with a tendency of bleeding. Due to limited published data and insufficient evidence of safety and efficacy of NM in CRRT, the acceptance of NM as a standard alternative anticoagulant has been hampered and the use of NM was initially limited to Japan. However, in 2005, NM was licensed in Korea when citrate anticoagulation is unavailable and the application of NM in CRRT has been increasing ever since (Ohtake et al., 1991;Nakae and Tajimi, 2003). It has been reported that NM attenuates ischemia reperfusion induced renal injury in mice, which is associated with the inhibition of apoptosis and reduction of nitric oxide overproduction mediated by NM (Na et al., 2016). However, the role of NM in the context of RIAKI is still unknown. Therefore, our study aimed to investigate the effect of NM acting rather than a simple anticoagulant on the RIAKI, and decipher the molecular mechanisms underlying the action of NM in RIAKI using integrated transcriptomic and metabolomic analysis.

Animals and drugs
Male Sprague Dawley (SD) rats (body weight: 226 g ± 13) were received from the animal experiment center of Southwest Medical University. The animals were housed in a facility with controlled humidity and temperature and 12 h dark/light cycle. Food and water were given ad libitum throughout the study. Glycerin was purchased from Soleibo Technology Co., Ltd. (Beijing, China). Nafamostat Mesylate was purchased from Durui Pharmaceutical Co., Ltd. (Jiangsu, China).

Experimental protocol
Thirty-six male SD rats were randomly divided into three groups: Control group (CN, n = 12), RIAKI group (RM, n = 12), and NM intervention group (NM, n = 12). The NM group was further divided into two subgroups: 24 h and 48 h observation groups (n = 6 in each group). The CN group received an intramuscular injection of normal saline (10 ml/kg); the RM group received an intramuscular injection of 50% (v/v) glycerol normal saline (10 ml/kg); and NM group received intraperitoneal injection of 1 mg/kg NM 30 min prior to the modeling. The CN group and the RM group also received an intraperitoneally injection of the same amount of normal saline.

Kidney function evaluation
Blood samples were collected without anticoagulants and centrifuged at 3,000 rpm for 10 min to harvest the serum. The serum blood urea nitrogen (BUN) and serum creatinine (SCr) were measured using an autoanalyzer (Siemens, Germany).

Histopathological evaluation
Kidney tissue was fixed in 10% (v/v) formalin, embedded in paraffin, and sectioned at 4 μm thickness. Following sectioning, H&E staining was performed. The extent of tubular vacuolation, dilatation, and necrosis were evaluated. Renal tubular damage was assessed using the following scoring system: 0 represents no tubular injury; 1 represents less than 25% tubular injury; 2 represents 25%-50% tubular injury; 3 represents 51%-75% tubular injury, and 4 represents more than 75% tubular injury. Five fields were randomly selected from each slide for assessments of renal tubular injury.
Terminal deoxynucleotide transferase dUTP nick end labeling staining Apoptosis in kidney sections was detected by terminal deoxynucleotide transferase dUTP nick end labeling (TUNEL) assay using the TUNEL apoptosis detection kit (Roche, Switzerland) according to the manufacturer's instructions. The number of TUNEL-positive cells and the total cells in 5 randomly selected fields of each slice were calculated as the apoptosis index (AI).

Measurement of mitochondrial membrane potential
Mitochondrial membrane potential (MMP) was detected by JC-1 dye (Haimen BiYunTian Biotechnology Research Institute, Nantong, China) according to the manufacturer's instructions. JC-1 aggregates in mitochondria in normal health cells, fluorescencing red, while in apoptotic cells, JC-1 accrues in the cytosol, as a green fluorescencing monomer. In brief, harvested cells from kidney tissue were incubated with 500 µl of JC-1 working solution in incubator for 20 min, washed, and resuspended with 500 µl of assay buffer, and analyzed by flowcytometry (cytoflex, Beckman).

Immunohistochemical staining
Kidney sections were deparaffinized, hydrated, and antigenretrieved. Endogenous peroxidase activity was quenched by incubating the specimens with a 3% (v/v) H 2 O 2 . The specimens were then blocked with a serum block solution (DAKO, Tokyo, Japan), followed by incubation with the Frontiers in Pharmacology frontiersin.org primary antibody (Rabbit monoclonal anti-GPX4; Cat. No. ab231174; Abcam) overnight at 4°C in a humidified box. After incubation with the secondary antibody (Goat Anti-Rabbit IgG ab6721, Abcam) for 30 min at room temperature, sections were incubated with DAB reagents for coloration (Abcam). Slides were examined using a Nanozoomer S60 digital slice scanning system. The positive areas of IHC staining (brownish yellow) were analyzed. Average optical density was used to assess the expression/signal intensity of GPX4.

Assessment of iron concentration and the GSH/GSSG ratio
The iron concentration of kidney tissue was measured by tissue iron assay kit (A039-2-1, Nanjing Jiancheng Bioengineering, China) according to the manufacturer's protocol. GSH/GSSG ratio of kidney tissue was determined using total glutathione/Oxidized glutathione assay kit (A061-2-1, Nanjing Jiancheng Bioengineering, China) according to the manufacturer's instruction.

RNA isolation and transcriptomic analysis
Total RNA was isolated using TRIzol reagent (Invitrogen, CA, United States). The RNA concentration and quality were evaluated using a Nanodrop 2000 instrument (Thermo Scientific, United States) and an Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, United States). The RNA-seq libraries were prepared using Illumina Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, CA, United States) according to the manufacturer's instructions. The transcriptome sequencing and analysis were conducted by OE Biotech Co., Ltd (Shanghai, China).
After removing the low-quality sequence reads and the adapters with fastp 1 , the clean reads were mapped to the rat genome using HISAT2 2 . FPKM 3 of each gene was calculated and the read counts of each gene were obtained by HTSeq-count 4 . PCA analysis was performed using R (v 3.2.0) to evaluate the biological duplication of samples. Differential expression analysis was performed using the DESeq2 5 . A p-value <0.05 and foldchange >2 or foldchange <0.5 was set as the threshold for significantly differential expression gene (DEGs). Hierarchical cluster analysis of DEGs was performed using R (v 3.2.0) to demonstrate the expression pattern of genes in different groups and samples. The radar map of the top 30 genes was drawn to show the expression of upregulated or downregulated DEGs using R packet radar.
Based on the hypergeometric distribution, Gene Ontology 6 and Kyoto Encyclopedia of Genes and Genomes (KEGG) 7 , pathway enrichment analysis of DEGs was performed to screen the significant enriched term using R (v 3.2.0), respectively. R (v 3.2.0) was used to draw the column diagram, the chord diagram and the bubble diagram of the significant enrichment term. Identification of transcription factors was conducted by matching DEGs with the rat transcription factors deposited in the Animal TFDB 3.0 database.

RT-qPCR analysis
The same RNA samples used for transcriptomic analysis were also used for RT-qPCR analysis. The yield of RNA was determined using a NanoDrop 2000 spectrophotometer (Thermo Scientific, United States), and the integrity was evaluated using agarose gel electrophoresis stained with ethidium bromide.
Quantification was performed with a two-step reaction process: reverse transcription (RT) and PCR. Each RT reaction consisted of 0.5 μg RNA, 2 μl of 5 × TransScript Allin-one SuperMix for qPCR, and 0.5 μl of gDNA remover, in a total volume of 10 μl. Reactions were performed in a GeneAmp ® PCR System 9700 (Applied Biosystems, United States) for 15 min at 42°C, and for 5 s at 85°C. The 10 μl RT reaction mix was then diluted ten times in nuclease-free water and held at −20°C. RT-qPCR was performed using LightCycler ® 480 II Real-time PCR Instrument (Roche, Swiss) with 10 μl PCR reaction mixture that included 1 μl of cDNA, 5 μl of 2 × PerfectStartTM Green qPCR SuperMix, 0.2 μl of forward primer, 0.2 μl of reverse primer, and 3.6 μl of nuclease-free H 2 O. Reactions were incubated in a 384well optical plate (Roche, Swiss) at 94°C for 30 s, followed by 45 cycles of 94°C for 5 s, and 60°C for 30 s. Each sample was run in triplicate for analysis. At the end of the PCR cycles, melting curve analysis was performed to validate the specific generation of the expected PCR product. Each sample was analyzed in three independent biological replicates, and the 2 −ΔΔCT method was used to determine the gene expression levels. The rat ACTB gene was used as an internal control. The primers were designed by OE Biotech Co., Ltd. (Shanghai, China), based on the mRNA sequences obtained from the NCBI database and synthesized by Tsingke Biotechnology Co., Ltd. (Beijing, China). Primers sequences are shown in Table 1.
The mass spectrometer was operated in positive and negative modes. The raw data was processed by Progenesis QI v2.4 (Waters, United States) following steps of alignment of the chromatograms, baseline filtering, normalization, integration, peak recognition, and retention time correction. The t-test and the score of variable importance in projection (VIP) were used to screen the differential metabolites between groups. The differential metabolites were mapped to the KEGG database for pathway enrichment analysis.

Combined metabolomic and transcriptomic analysis
Pearson correlation coefficients and p values were used to screen metabolites and related genes within the conjoint analysis of metabolomics and transcriptomics. To further understand the relationship between DEGs and differential metabolites, the DEGs and metabolites were mapped to their associated KEGG pathways.

Statistical analysis
Statistical analysis was performed using IBM SPSS 25.0 statistical software (IBM Corp.). The data were expressed as mean ± standard deviation comparison between groups was performed by one-way ANOVA, and pairwise comparison between groups was performed by S-N-K test. A Kruskal-Wallis test was performed when data were heteroscedastic and a result of p < 0.05 was considered statistically significant. For the transcriptomic analysis, the DEGs were detected with a false discovery rate-adjusted p-value <0.05 and foldchange >2 or foldchange <0.5. For the metabolomic analysis, the threshold for the t-test was p < 0.05, and the threshold for VIP score was set at ≥ 1.

Results
The effect of nafamostat mesylate on renal function in rhabdomyolysis-induced acute kidney injury There was a significant increase in serum BUN and SCr in the RM group at both 24 h and 48 h compared to the CN group (p < 0.01). The levels of serum BUN and SCr decreased significantly in the NM group compared to the RM group at 24 h (p < 0.01), and there was a significant reduction in serum BUN and SCr following 48 h intervention of NM (p < 0.01) (Figure 1).

Effect of nafamostat mesylate on inflammatory cytokines and proenkephalin in rhabdomyolysis-induced acute kidney injury
There was a significant increase in the serum levels of IL-6, TNF-α and PENKID in the RM group compared to the CN group at both 24 h and 48 h after the glycerol challenge (p < 0.01). The level of serum IL-6 showed the trend of dropping at 48 h, while the level of serum TNF-α and PENKID continued to increase overtime. Nafamostat mesylate intervention significantly reduced the serum levels of IL-6, TNF-α and PENKID compared to the RM group (p < 0.01) (Figure 2).

Effect of nafamostat mesylate on renal histopathology in rhabdomyolysisinduced acute kidney injury
In H&E staining, a higher degree of vacuolization and tubular dilatation, and a large number of necrotic cell and myoglobin casts were noted in the RM group, whereas this phenomenon was mitigated with NM intervention. The renal tubular injury scores were significantly increased in the RM group compared with the CN group at both 24 h and 48 h time points (p < 0.01), while there was a significant decrease in renal tubular injury scores in the NM group compared with the RM group at 24 h (p < 0.05) and 48 h (p < 0.01) (Figure 3).

Effect of nafamostat mesylate on renal tissue apoptosis in rhabdomyolysisinduced acute kidney injury
There were more TUNEL-positive cells in the RM group compared to the CN group at both 24 h and 48 h, while fewer TUNEL-positive cells were seen in the NM group compared with the RM group ( Figure 4A). The AI staining showed that apoptosis was significantly increased in the RM group compared to the CN group at both 24 h and 48 h time points Frontiers in Pharmacology frontiersin.org 06 (p < 0.01), whereas it was mitigated significantly with NM intervention at both time points (p < 0.01) ( Figure 4B).

Effect of nafamostat mesylate on renal tissue ferroptosis in rhabdomyolysisinduced acute kidney injury
The iron accumulations increased significantly in the RM group at 24 h and 48 h compared to CN group at the corresponding time points (p < 0.05), and decreased after NM treatment at 24 h compared to RM group (p < 0.05). There is no significant difference in kidney tissue iron contents between NM and RM group at 48 h. The GSH/GSSG ratio showed a remarkable decrease in RM group compared to CN group at both time points (p < 0.01), while there was a significant increase in GSH/GSSG ratio after RM treatment at 24 h and 48 h time points compared to RM group (p < 0.05) (Figures 5A,B).
The expression of GPX4 was lower in the RM group compared to the CN group at both 24 h and 48 h (p < 0.01), while the abundance of GPX4 increased significantly after NM treatment (NM vs. RM, p < 0.01) (Figures 5C,D). The RM group showed an enhanced loss of MMP compared with the CN group, while NM treatment showed an improved loss of MMP compared with the RM group (p < 0.05)  Frontiers in Pharmacology frontiersin.org ( Figures 5E,F). This suggests that the NM attenuates the loss in MMP in response to RIAKI.

Differential gene expression profiles in rhabdomyolysis-induced acute kidney injury
The screening criteria for DEGs were a false discovery rateadjusted p ≤ 0.05 and log2 fold change >1. The transcriptome analysis identified 3,891 DEGs in the RM group compared with the CN group, with 2,132 genes upregulated and 1,759 genes downregulated, and 294 DEGs in the NM group compared with the RM group, of which 182 genes were upregulated and 112 genes were downregulated ( Figure 6). To gain further insight into related signaling pathways, KEGG pathway analysis was performed and revealed that 192 signaling pathways were enriched in the NM group compared with the RM group, of which 147 were up-regulated and 103 were down-regulated. The upregulated pathways were mainly implicated in glutathione metabolism, steroid biology synthesis, amino acid (tyrosine, glycine, serine, threonine, pyruvate, tryptophan, cysteine, and methionine) metabolism, arachidonic acid metabolism, metabolism of cytochrome P450-related substances, and peroxisome. The down-regulated pathways were classified in IL-17 signaling pathway, cytokine-cytokine receptor interaction, PPAR signaling pathway, chemokine signaling pathway, tumor necrosis factor signaling pathway, and bile secretion ( Figures 7A,B). Five hundred and five DEGs (291 up and 214 down) were annotated with the transcription factor database. Compared with the RM group, the most of upregulated DEGs in NM group belonged to Glutathione metabolism, whereas the most of downregulated DEGs belonged to the transcription factor Cytokine-cytokine receptor interaction ( Figure 7C). Notably, the glutathione and IL-17 signaling pathways were assigned the most DEGs.

Verification of the key DEGs profiles in rhabdomyolysis-induced acute kidney injury
To further verify the credibility of the expression profiles of the key DEGs obtained from the RNA-seq data, 7 DEGs were selected based on the most enriched signaling pathways and a large fold change. The RT-PCR results showed that the expression trends of these DEGs were highly consistent with the transcriptome data. The nafamostat mesylate mediated glutathione metabolic pathway related genes Anpep, Gclc, Ggt1 were upregulated, while Mgst2 was significantly downregulated (p < 0.01). The inflammatory pathway related genes Cxcl13 and Rgn were significantly upregulated (p < 0.01) and the ferroptosis pathway related gene Akr1c1 was significantly upregulated (p < 0.01) (Table 2; Figure 8). Collectively, these results suggest that various biological pathways are involved in renal protection with NM in RIAKI.

Correlation of metabolic changes and gene expression in rhabdomyolysisinduced acute kidney injury
The metabolomic analysis identified a total of 559 differential metabolites in the RM group compared with Frontiers in Pharmacology frontiersin.org the CN group, and 262 differential metabolites with differential changes under NM intervention compared with RM. The KEGG enrichment analysis of differential metabolites showed the enriched differential metabolites were mainly involved in amino acid metabolism (lysine degradation, valine, leucine and isoleucine biosynthesis, arginine and Frontiers in Pharmacology frontiersin.org 09 proline metabolism, D-arginine and D-ornithine metabolism), choline metabolism in cancer, mTOR signaling pathway, glycerophospholipid metabolism, necrosis, apoptosis, synthesis and secretion of cortisol, vitamin digestion and absorption, and pyrimidine metabolism (Figure 9). Meanwhile, the integrated transcriptomics and metabolomics analyses revealed the association network of key DEGs and corresponding metabolites, among them, Gclc interacts with 16 differential metabolites, of which 8 metabolites were positively correlated with Gclc, while 8 metabolites were negatively correlated with Gclc (Table 3; Figure 10). This data suggests that the accumulation of Gclc may contribute to renal protection with NM in RIAKI.

Discussion
Our study demonstrated that NM exerts a protective effect against AKI in the context of RM. The further integrated transcriptomic and metabolomic analysis revealed that the DEGs regulated by NM were mainly implicated in the glutathione metabolism pathway, inflammatory response, and ferroptosis related pathways. More importantly, these pathways related genes Anpep, Gclc, Ggt1, Mgst2, Cxcl1, Cxcl3, Rgn, Retn, Akr1c1, and Akr1c2 were highly regulated by NM. Based on our results, it can be speculated that the renal protective effect of NM on RIAKI is mainly produced through the regulation of glutathione metabolism, inflammatory response, and ferroptosis. Targeting these key pathways related genes may emerge as a potential molecular therapy for RIAKI.
Rhabdomyolysis is a common cause of AKI, and the mechanisms involved in RIAKI are complex. Therefore, further understanding of the potential molecular mechanism of RIAKI and the identification of new targets will facilitate the exploring of novel therapies for RIAKI. In the present study, we developed a glycerol induced experimental model of RIAKI in rats using well established methods (Al Asmari et al., 2017). Our model successfully replicated the pathogenesis of RIAKI characterized by the deterioration of kidney function, inflammation, cell apoptosis, and histopathological changes. We observed a significant increase in serum SCr and BUN over time after modeling, which indicated the progressive deterioration of kidney function. There was a significant increase in inflammatory cytokines and cell apoptosis in the modeling group. Our histological data also showed tubular pathological damages with a large amount of necrotic cells and myoglobin. Myoglobin is known to play an important role in the development of RIAKI (Gburek et al., 2003). Myoglobin is an iron and oxygen binding protein found generally in the muscle tissue of vertebrates and in almost all mammals. It has a higher binding affinity for oxygen than hemoglobin (David, 2000). Increased evidence has suggested that released myoglobin in circulation can activate the reninangiotensin-aldosterone system, and trigger the release of vasopressin and nitric oxide deficiency, resulting in intrarenal vasoconstriction, tubular obstruction and direct tubule injury, and inflammation (Holt and Moore, 2001;Belliere et al., 2015). During RM, excessive amounts of myoglobin are released into circulation from muscle cells, and the free myoglobin is filtered by the glomerular filtration barrier and endocytosed by renal Frontiers in Pharmacology frontiersin.org tubular cells via the epithelial endocytic receptors megalin and cubilin (Khan, 2009). Ferrous myoglobin inside tubular cells is oxidized to a ferric form, resulting in the production of a hydroxyl radical, followed by the transformation of ferric myoglobin into ferryl myoglobin through redox cycling, yielding of a radical species. These radical species lead to lipid peroxidation and excessive oxidative stress (Petejova and Martinek, 2014). The clinical treatment for RIAKI is limited. The standard medical interventions are usually the management of underlying disease and initiation of conservative treatments, including massive hydration, urine alkalization and forced diuresis (Krouzecky et al., 2003). Currently, the core treatment is renal replacement therapy (RRT). It has been reported about 8%-65% of patients who develop RIAKI are in need of RRT (de Meijer et al., 2003). The possibility of myoglobin removal using either intermittent hemodialysis or CRRT in the case of RIAKI has been demonstrated in many studies (Amyot et al., 1999;Heyne et al., 2012). It has been shown in several retrospective studies that NM is safe and effective for patients with high bleeding risk in CRRT (Maruyama et al., 2011;Baek et al., 2012;Hwang et al., 2013), and does not affect patient mortality (Baek et al., 2012;Lee et al., 2014). Choi et al. Choi et al. (2015) reported that patients who received NM treatment had better survival rates at 30 and 90 days after CRRT initiation than when compared to the non-NM treatment group, however, the difference was not statistically Frontiers in Pharmacology frontiersin.org significant. So far, there is very limited data regarding the protective effect of NM on RIAKI besides acting as an anticoagulant in CRRT. Our study found, for the first time, that NM intervention offers a beneficial effect in RIAKI in a rat model. As is known there is excessive oxidation during RM induced renal injury (Moore et al., 1998;Giannoglou et al., 2007). Glutathione (GSH) metabolism is one of the most important antioxidant systems against oxidative stress, which plays an important role in regulating oxidative stress in many diseases, such as diabetes, stroke, and neurodegenerative disease etc. (Wu et al., 2004). The downregulation or deficiency of GSH has been associated with excessive oxidative stress in renal ischemic repercussion disease (Shang et al., 2016). Glutathione protects cellular components against lipid peroxidation, and DNA and protein damage through neutralizing hydrogen peroxide and reactive oxygen species (ROS) (Hanschmann et al., 2013) Cellular GSH is mostly present in the cytosol and can be oxidized by ROS into glutathione disulfide (GSSG). Therefore, the ratio of GSH/GSSG can indicate the cellular redox state. Glutamate-cysteine ligase (GCL) is the rate limiting enzyme for GSH synthesis. It consists of a catalytic subunit Gclc and a modifier subunit Gclm (Giordano et al., 2007). The catalytic Gclc subunit is essential for the enzymatic activity, while the modulatory Gclm subunit enhances Gclc catalytic efficiency. It has been reported that a lack of GCLM results in a reduction of GSH and consequent vulnerability to oxidative stress. The expression of GCL is induced by oxidative stress and this process can be regulated by some pivotal transcription factors, such as activator protein-1, nuclear factor-erythroid 2 related factor 2, and nuclear factor κB (NFκB) etc. (Lu, 2009). The amount of GCL subunits, the availability of substrates, and the extent of feedback inhibition of GCL by GSH control the overall rate of GSH synthesis. Botta et al. Botta et al. (2004) reported that overexpression of Gclc and GCLM in mouse liver hepatoma cells exhibits increased resistance to TNF-induced mitochondria injury and apoptosis. Tran's study Tran et al. (2004) showed that adenoviral overexpression of the GSH catalytic subunit Gclc increases intracellular GSH levels and protects beta cells against oxidative stress. Therefore, the strategy favors the biosynthesis of GSH and can be an attractive antioxidant therapy. Boutaud et al. Boutaud et al. (2010) reported that acetaminophen attenuates RM induced renal failure by reducing radical species generated from redox cycling between ferric and ferryl myoglobin. Chander et al. Chander et al. (2003) showed that a natural antioxidant, catechin, can reduce the toxicity of circulating myoglobin in the kidney tissues, and protect the kidneys against myoglobinuria acute renal failure in rats. In line with these studies, our results found that NM protects the kidneys against RIAKI, which may be associated with the upregulation of Gclc through activating the glutathione metabolism pathway. This hypothesis was supported by our observation that NM treatment increased GSH/GSSG ratio significantly compared to RM group. In the meantime, we observed an increased expression of Ggt1 after NM intervention, which is consistent with a report that showed the protective action of GSH requires GSH hydrolysis catalyzed by Ggt1 (Paolicchi et al., 2003). Our conjoint analysis of the transcriptome and metabolome also showed that Gclc interacted with 16 key differential metabolites, which suggests Gclc may be the key gene implicated in the renal protective effect of NM in RIAKI. This provides a mechanistic basis for antioxidant therapy in RIAKI. The physiological function of Mgst2 as a mediator of oxidative stress and DNA damage in nonimmune cells has been well studied. Dvash et al. (2015)'s study has shown that Mgst2 deficiency attenuates ER stresstriggered oxidative stress, DNA damage and apoptosis, and mouse morbidity in acute kidney injury. Consistent with this study, our data showed that NM treatment decreases the transcription level of Mgst2, which may contribute to renal protection. It is known that inflammation plays a pivotal role in the pathogenesis of RIAKI and that the leukocytes, particularly macrophages, are involved in the development of RIAKI (Holt and Moore, 2001). Macrophage infiltration has been reported in RM patient's kidney biopsy specimens and glycerol-induced AKI in rat models (Sinniah and Lye, 2000;Homsi et al., 2006). Myoglobin can stimulate proximal tubular cells to secrete macrophages and many molecules, such as high mobility group box proteins, uric acid, and microRNAs. These molecules are released from damaged muscle cells and activate the immune system, which lead to the activation of dendritic cells, Toll-like receptors, T-lymphocytes, macrophages, and NFκβ, and result in the increased generation of pro-inflammatory cytokines, such as IL-1 and TNF. Increased expression of adhesive molecules (ICAM-1, VCAM-1 etc.) also enhances the inflammatory response. Myoglobin-derived heme is another important molecular that influences the inflammatory reaction in the tubular endothelium and epithelium cells (Jang and Rabb, 2009;Panizo et al., 2015). Therefore, anti-inflammatory treatment can be a potential therapeutic strategy in RIAKI. Macrophage depletion has been reported to protect the kidneys against ischemia-reperfusion induced AKI (Ko et al., 2008). A phosphodiesterase inhibitor, pentoxifylline, was shown to decrease inflammation through the inhibition of TNF-α and leukotriene production in RIAKI (Plotnikov et al., 2009). In agreement with the above studies, our results showed that NM offers kidney protection through reducing the proinflammatory cytokines IL-6 and TNF-α in RIAKI. Our data showed that NM also upregulated the anti-inflammatory gene Rgn, which is consistent with the report that Rgn may act as a suppressor of inflammation in macrophage-infiltrated adipocyte tissue (Murata et al., 2020). Surprisingly, we found that NM induced the expression of the pro-inflammatory chemokine Cxcl13. However, the underlying cause is still not clear. Nevertheless, the effect of NM induced upregulation of Cxcl13 was overwhelmed by the NM mediated renoprotection. In the setting of RM, lipid peroxidation and oxidative stress can stimulate the release of cytochrome c from the mitochondria which lead to the activation of caspase 1 and 3, resulting in apoptosis (Panizo et al., 2015). However, other types of cell death may also be involved in RIAKI. It has been reported that there was no change in creatinine levels in mice that received specific apoptosis inhibitors (Homsi et al., 2006). Guerrero-Hue et al. (2019) demonstrated another type of nonapoptotic cell death. Ferroptosis plays a pivotal role in kidney injury associated with RM in vitro myoglobin stimulated renal tubular cells and in the mouse model of RIAKI, and Melania et al. found that ferroptosis is consistent with the increased lipid peroxidation and iron accumulation in the context of RM. More importantly, they found that the inhibition of ferroptosis attenuated kidney dysfunction, whereas no effects were observed with the inhibition of apoptosis or necroptosis, which implicates a FIGURE 8 mRNA expression of key DEGs in RIAKI (*p < 0.05, **p < 0.01).
Frontiers in Pharmacology frontiersin.org prominent contribution of ferroptosis rather than apoptosis to the RIAKI. Many other studies have also reported that ferroptosis contributes to kidney injury caused by cisplatin, folic acid, and ischemia-reperfusion injury in experimental models (Linkermann et al., 2014a;Linkermann et al., 2014b;Martin-Sanchez et al., 2017). Our study investigated the ferroptosis related protein GPX4 in the context of RIAKI. As is known, GPX4 is a key GSH-utilizing enzyme which counteracts lipoxygenase activity and phospholipid/ cardiolipin oxidation (Conrad and Friedmann Angeli, 2015). Inhibition of GPX4 causes the accumulation of iron-dependent lipid peroxidation, which ultimately leads to ferroptotic cell death (Chen et al., 2020). Friedmann Angeli et al. (2014) reported that GPX4 knockout in tubular cells triggers ferroptosis and substantial cell death in acute kidney injury. In line with these studies, our data showed the repression of GPX4 and increased accumulation of iron in RIAKI, while the expression of GPX4 increased significantly after NM treatment accompanied with decreased iron levels, which supports the idea that ferroptosis may also be involved in RIAKI and can be regulated by NM. Meanwhile, we observed that NM intervention mitigates apoptosis induced by RIAKI. Considering the TUNEL staining used in our study is not specific for apoptosis, we also investigated the effect of NM on MMP, and found NM treatment attenuates the loss in MMP in response to RIAKI. The aldo-keto reductases (AKRs) gene family has been reported to inhibit ferroptosis through potentiating the metabolization of lipid peroxides by catalyzing aldehydes and FIGURE 9 (A) Volcano plot of differential metabolites in RIAKI. X axis and Y axis represent log 2 (fold change) and p-value, respectively. Red and green dots represent upregulated and downregulated metabolites, respectively. (B) KEGG pathway enrichment analyses of differential metabolites. (B) The top 20 are presented. The enrichment factor is the ratio of the number of DEGs to the number of total genes in a particular pathway. The dot size represents the number of differential metabolites, while colors correspond to the range of p-values.
Frontiers in Pharmacology frontiersin.org ketones to alcohols (Dixon et al., 2012). Gagliardi et al. (2019) study showed that AKR1C1/3 degrades the 12/15-LOXgenerated lipid peroxides, thus resulting in ferroptotic cell death resistance. It has been reported that AKR1C1 protects cancer cells from ferroptosis and is considered to be the ferroptosis-protective gene (Wohlhieter et al., 2020). These studies support the notion that AKR1C1 plays a role in ferroptosis. Our data also showed that AKR Family one Member C1 (Akr1c1) and Akr1c2 were both highly regulated by NM, which suggests that the inhibition of ferroptosis contributes to the NM mediated renal protection in RIAKI. We appreciate the limitations of our study. First, we treated the RIAKI rats with NM at a dose of 1 mg/kg and did not examine the dose-response relationship of the effect of NM on RIAKI. Second, we administered NM prior to the induction of RIAKI, which is not amenable in clinical practice. Thus, the therapeutic effect of NM on RIAKI remains to be elucidated in further studies. Another limitation is that the sample size of our study was relatively small, which may cause a certain degree of sampling bias. Finally, our study did not define the molecular functions of the key DEGs and metabolites identified and their possible mechanisms by which RIAKI is attenuated with treatment of NM. Rather, our work provides the transcriptomic and metabolomic profiling and framework from which the key molecules can be prioritized for future mechanistic studies.

Conclusion
Our study comprehensively analyzed the changes in the abundance of genes and metabolites regulated by NM during RIAKI. The results revealed that NM broadly activates glutathione metabolism, inflammatory, and ferroptosis associated signaling pathway, which leads to accumulation of antioxidant, anti-inflammatory, and anti-ferroptosis related metabolites and genes. Our findings provided insight into the regulatory mechanism of NM mediated renal protection, which facilitates the exploration and translation of therapeutic or diagnostic properties of these candidate genes to the clinic. However, the functional characterization of the identified DEGs in this study still needs to be further investigated.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/bioproject/PRJNA851578.

Ethics statement
The animal study was reviewed and approved by experimental animal ethics committee of Southwest Medical University.

FIGURE 10
Conjoint analysis of transcriptome and metabolome in RIAKI. Green and red circles represent corresponding DEGs and metabolites, respectively. Green and red lines indicate positive and negative correlations between DEGs and metabolites, respectively.
Frontiers in Pharmacology frontiersin.org