Integrated Analysis of m6A Methylome in Cisplatin-Induced Acute Kidney Injury and Berberine Alleviation in Mouse

Background N6-methyladenosine (m6A) is the most abundant modification known in mRNAs. It participates in a variety of physiological and pathological processes, such as metabolism, inflammation, and apoptosis. Aims To explore the mechanism of m6A in cisplatin-induced acute kidney injury (AKI) and berberine alleviation in mouse. Methods This study investigated the N6-methyladenosine (m6A) methylome of kidneys from three mouse groups: C57 mice (controls), those with CI-AKI (injury group, IG), and those pretreated with berberine (treatment group, TG). Methylated RNA Immunoprecipitation Next Generation Sequencing (MeRIP-seq) and RNA-seq were performed to identify the differences between the injury group and the control group (IvC) and between the treatment group and the injury group (TvI). Western blotting was performed to identify the protein levels of candidate genes. Results In IvC, differentially methylated genes (DMGs) were enriched in metabolic processes and cell death. In TvI, DMGs were enriched in tissue development. Several genes involved in important pathways related to CI-AKI showed opposite methylation and expression trends in the IvC and TvI comparisons. Conclusion m6A plays an important role in cisplatin induced AKI and berberine may alleviate this process.


INTRODUCTION
Cisplatin is an anticancer drug widely used for the treatment of various solid tumors, but it is also known for its nephrotoxicity. It exerts its function by interacting with and disrupting DNA and mitochondrial function. During drug metabolism, cisplatin accumulates in renal tubular cells, causing cell death and resulting in acute kidney injury (AKI) (Lebwohl and Canetta, 1998;Cummings and Schnellmann, 2002;Ozkok and Edelstein, 2014;George et al., 2018;Holditch et al., 2019;Yimit et al., 2019). Recent studies have revealed that apoptosis, necrosis, inflammation, and other mechanisms play significant roles in cisplatin-induced AKI (Sahu et al., 2015;Zuk and Bonventre, 2016;Humanes et al., 2017;Long et al., 2017).
m6A has been implicated in all aspects of posttranscriptional RNA metabolism, such as regulating reversible modifications, alternative splicing, stability, and translation. It prefers to modify sequences identified as RRACH, where R is adenine or guanine, A is the m6A site, and H is adenine, cytosine, or uracil. In addition, m6A modifications exhibit enrichment in the 3' UTR near mRNA stop codons and within long internal exons. The writing of m6A is accomplished via a complex composed of methyltransferase-like 3 (METTL3), methyltransferase-like 14 (METTL14), and Wilms tumor 1-associated protein (Ping et al., 2014;Wang et al., 2014;Spitale et al., 2015). Two m6A demethylases, fat mass and obesity-associated protein (FTO) and AlkB homolog 5, have been discovered as "erasers" (Jia et al., 2011). While proteins containing the YT521-B homology (YTH) domain, such as YTHDC1, YTHDF1, YTHDF2, and YTHDF3, directly bind m6A sites and act as readers of the m6A signal (Zheng et al., 2013).
Several studies have focused on the correlation between m6A methylation and kidney injury. Wang J. et al. (2019) reported that METTL3 overexpression plays a protective role against colistin-induced oxidative stress and apoptosis in renal tubular epithelial cells in mice. Alteration of m6A regulators was associated with pathologic stage in patients with clear cell renal cell carcinoma (Zhou J. C. et al., 2019). Another study reported that cisplatin treatment reduced FTO expression and increased m6A levels in vivo and in vitro. They also found that inhibiting FTO by meclofenamic acid aggravated renal damage and increased apoptosis in cisplatin-treated kidneys (Zhou P. et al., 2019). METTL14 is elevated in people with AKI .
However, few studies have investigated the m6A methylome in cisplatin-induced AKI and the potential mechanism of berberine alleviation.
In this study, we found that berberine significantly alleviated cisplatin-induced AKI in a reliable mouse model. To further investigate the role of m6A in this process, MeRIP-seq was used to establish the first known transcriptome-wide m6A methylome profiles of kidneys from normal, CI-AKI, and berberine-pretreated mice. RNA-seq was performed to detect differentially expressed genes (DEGs) among the groups. Based on our results, we speculate that berberine may alleviate CI-AKI by regulating m6A methylation.

Animals and Tissue Collection
Male C57BL/6 mice (aged 8 weeks) were randomly assigned to the control group, injury group (IG), and treatment group (TG), with four mice per group. All mice were housed under a 12 h light/dark schedule with free access to food and water. Control and IG mice were subjected to daily intraperitoneal (i.p) injections with vehicle (normal saline), while TG mice were injected daily with berberine (Sigma, St. Louis, MO, United States, 20 mg/kg) (Ruan et al., 2017). After 14 days of drug treatment, the IG and TG were injected intravenously with cisplatin (20 mg/kg), while controls were injected intravenously with a vehicle. After cisplatin injection, berberine pretreatment was stopped and all mice were housed as usual. All mice in the three groups were sacrificed at 72 h postinjection by cervical dislocation after CO 2 -induced narcosis (Zhang et al., 2016;Dutta et al., 2017). Immediately afterward, the kidneys were collected.

Serum Levels of Creatinine and Blood Levels of Urea Nitrogen
Before the mice were sacrificed, blood was drawn from their tail veins. Serum samples were collected. Serum levels of creatinine (Scr) and urea nitrogen (BUN) were analyzed using a standard spectrophotometric assay (Roche Diagnostics GmbH, Mannheim, Germany).

Histopathology Analyses
Renal tissue harvested from animals was washed with 0.9% saline, fixed in 10% neutral buffered formalin, and then embedded in 10% paraffin. Sections (5 µm thick) were stained with hematoxylin eosin for further microscopic analyses. A tubular injury score was calculated (Leemans et al., 2005). The percentage of damaged tubules in the corticomedullary junction was estimated by a nephropathologist using a 5-point scale according to the following criteria: tubular dilation, cast deposition, brush border loss, and necrosis in eight randomly chosen, nonoverlapping fields (×400 magnification). Lesions were graded on a scale from 0 to 5: 0 = normal; 1 = mild, involvement of less than 10% of the cortex; 2 = moderate, involvement of 10-25% of the cortex; 3 = severe, involvement of 25-50% of the cortex; 4 = very severe, involvement of 50-75% of cortex; 5 = extensive damage, involvement of more than 75% of the cortex.

RNA MeRIP-seq and Data Analyses
In accordance with the manufacturer's instructions, TRIzol reagent (Invitrogen Corporation, CA, United States) was used to extract total RNA from kidney tissue. Ribosomal RNA was removed from total RNA with a Ribo-Zero rRNA Removal Kit (Illumina, Inc., CA, United States). Then, fragmentation buffer (Illumina, Inc.) was used to split the RNA into fragments of approximately 100 nucleotides in length. MeRIP-Seq was performed by Cloudseq Biotech Inc. (Shanghai, China). In brief, total RNA was extracted from kidney tissue using TRIzol Reagent (Life Technologies CA, United States). RNA was tested for quality via NanoDrop (Thermo Fisher Scientific, MA, United States). RNA integrity was assessed using a denaturing agarose gel. mRNA was isolated from total RNA using a Seq-Star TM poly(A) mRNA Isolation Kit (Arraystar, MD, United States). The GenSeq TM m6A RNA IP Kit (GenSeq Inc., China) was used to perform m6A RNA immunoprecipitation. A NEBNext R Ultra II Directional RNA Library Prep Kit (New England Biolabs, Inc., MA, United States) was used to construct both the input samples without immunoprecipitation and the m6A IP samples. All samples were subjected to 150 bp paired-end sequencing on an Illumina HiSeq instrument (Illumina, Inc.).
Frontiers in Genetics | www.frontiersin.org P < 0.05 were identified using the diffReps differential analysis package (Shen et al., 2013;Wang Y. et al., 2019;Wang et al., 2020). The peaks identified by MACS and diffReps overlapping with exons of mRNA were selected for further analyses. Gene ontology (GO) and pathway enrichment analyses were performed on the differentially methylated proteins for using the GO 1 and Kyoto Encyclopedia of Genes and Genomes (KEGG) 2 databases.

RNA-Seq and Data Analyses
Total RNA was extracted from kidney samples using TRIzol reagent (Life Technologies) according to the manufacturer's protocol. Denaturing agarose gel electrophoresis was used to evaluate the integrity of total RNA. A Seq-Star TM poly(A) mRNA Isolation Kit (Arraystar, MD, United States) was utilized to purify mRNA from total RNA after measuring 1 www.geneontology.org 2 www.genome.jp/kegg the quantity and quality on a NanoDrop ND-2,000. Then, a BGISEQ-500 platform was used to subject fragmented mRNA to 50 bp single-end sequencing. Adapter and lowquality reads were trimmed using SOAPnuke (Chen et al., 2018), and those trimmed reads were aligned to the reference genome using bowtie2 (Langmead and Salzberg, 2012). Finally, cuffdiff was use to analyze differential expressed genes (DEGs) (Trapnell et al., 2012).

Western Blotting
Mouse tissues were lysed using a protein lysis buffer containing 20 mM Tris (pH 7.4), 150 mM NaCl, 1 mM EDTA, 1 mM EGTA, 1% Triton X-100, 25 mM sodium pyrophosphate, and 2 mM sodium orthovanadate aprotinin. All denatured proteins were separated on an SDS-PAGE gel and then transferred to polyvinylidene difluoride membranes (Roche, Netley, NJ, United States). The membranes were blocked with 5% skimmed milk in Tris-buffered saline and then were incubated Frontiers in Genetics | www.frontiersin.org with 1:500 dilutions of primary antibodies as follows: anti-FGA (Abcam, Cambridge, MA, United States), anti-Slc12a1 (Abcam, Cambridge, MA, United States), and anti-Havcr1 (Abcam, Cambridge, MA, United States). Then, the samples were incubated with a horseradish peroxidase-conjugated antirabbit secondary antibody (Jackson ImmunoResearch, PA, United States). The bands were visualized using an ECL Western Blotting Kit (Biovision, Milpitas, CA, United States) and were quantified by Quantity One software (Bio-Rad, Hercules, CA, United States).

Statistical Analyses
Data are expressed as mean ± standard deviation (SD). Student's t-test was used to compare two groups. ANOVA with Tukey's post-test was used to assess the statistical significance betweengroup means for comparisons between multiple groups. The differences were considered statistically significant at P < 0.05.

Establishment of a Reliable CI-AKI Mouse Model
We found that 20 mg/kg cisplatin injection resulted in an approximately 18-20-fold increase in Scr and BUN relative to the control group, but significantly lower levels in TG than in IG (Figures 1A,B). Hematoxylin-eosin staining indicated that cisplatin induced remarkable renal structure damage in IG tissue, including extensive tubular vacuolization, tubular epithelial cell exfoliation, and thickening of glomerular basement membrane (P < 0.001), while these changes were notably alleviated in TG mice (P < 0.001) (Figures 1C,D). This indicated successful and reliable establishment of a CI-AKI mouse model and relief of the injury through berberine pretreatment.

General Features of m6A Methylation
MeRIP-seq analyses of mRNA derived from kidneys revealed 13,284 m6A peaks within 7,942 coding gene transcripts in controls; the values were 11,846 within 7,422 mRNAs in IG and 10,106 within 6,481 mRNAs in TG. Overall, 7,337 peaks overlapped among the three groups (Supplementary Figure 1). Pairwise comparisons showed that 9,008 peaks (more than 55.8% of all peaks) overlapped in IG versus control (IvC) comparisons and 8,537 peaks (more than 63.6% of all peaks) overlapped in TG versus IG (TvI) comparisons (Figure 2A). Approximately 40% of all peaks were non-overlapping, suggesting differences among those groups.
To determine whether the m6A peaks contained the RRACH conserved sequence motif, we detected all samples in three groups. One thousand peaks with the highest scores (−10 * log10, p-value) were analyzed using DREME software. The results showed that the same sequence motif (RRACH) was necessary for m6A methylation in kidney mRNAs in each group (Figure 2B), consistent with other studies (Dominissini et al., 2012;Meyer et al., 2012;Luo et al., 2019).
Further analyses demonstrated that 16.59, 21.21, and 22.29% of the genes with m6A-methylation sites in controls, IG, and TG contained one m6A peak ( Figure 2C). More than 75% of the genes contained two or more peaks ( Figure 2C). This result is different from other reports on mouse liver (Luo et al., 2019) and brain (Dominissini et al., 2012), suggesting that the kidney is unique.
To understand the preferential locations of m6A peaks, all peaks were categorized into five transcript segments: the 5'  UTR, start codon segment, coding sequence (CDS), stop codon segment, and 3' UTR. We found that m6A was mostly enriched in the CDS, and there was some enrichment near stop codons ( Figure 2D). These results are also similar previous studies (Dominissini et al., 2012;Meyer et al., 2012). Furthermore, m6A peaks in all groups had the highest density in stop codons ( Figure 2E). Comparisons between different groups showed no significant difference in the volume of m6A peaks, which indicates that although many variants in m6A methylation were detected in different segments, the total proportion of m6A peaks did not significantly change ( Figure 2F).
To understand the DMMS distribution profiles in those two groups, we mapped them to chromosomes. Chromosomes 4, 2, 7, and 11 were the top four chromosomes, harboring more than 200 DMMSs (Supplementary Figure 2A). However, when the number of DMMSs was normalized by the length of chromosomes, the top four with the highest relative densities were 11, 19, 7, and 4 ( Figure 3A). Meanwhile, most DMMSs were within a CDS (Figure 3B). Further analyses showed that at sites with increased methylation, DMMSs within the 5' UTR had the highest fold change, while among the sites with decreased methylation, DMMSs within the 3' UTR had the highest fold change (Figure 3C), demonstrating the location preferences of methylation and demethylation within the genome.
In TvI comparisons, 526 DMMSs were identified within 420 genes, of which 66.73% (351/526) were sites with increased methylation. The top 20 m6A sites with the most increased and decreased methylation are shown in Table 2.
Chromosomes 3, 1, 9, and 4 were the top four chromosomes harboring the most DMMSs (Supplementary Figure 2B). The top four with the highest relative densities were 9, 19, 18, and 3 ( Figure 3D). Most DMMSs were within a CDS (Figure 3E). DMMSs within the 3' UTR had the highest fold change among sites with increased methylation, while those near the start   codon had the highest fold change among sites with decreased methylation (Figure 3F). We further analyzed DMMSs with contrary methylation trends in IvC and TvI. In total, 94 DMMSs showed opposite trends. Of these, 42.55% (40/94) exhibited significantly increased methylation in IvC but decreased methylation in TvI.

Differentially Methylated RNAs Are Involved in Important Biological Pathways
To explore the role of m6A in different groups, GO enrichment analyses and KEGG pathway analyses were used to analyze all protein coding genes containing DMMSs.
In IvC (Figure 4), for the BP category, genes with increased methylation of m6A sites were significantly (p < 0.05) enriched in the regulation of metabolic processes and cell death processes, such as macromolecule metabolic, cellular metabolic, and apoptotic processes ( Figure 4A); genes with decreased methylation were highly enriched in metabolic processes, oxidation-reduction processes, transport, transmembrane transport, and others ( Figure 4B). Regarding pathways, the former were significantly involved in apoptosis-associated pathways (e.g., TNF, MAPK, P53 signaling pathways) while the latter were involved in propanoate metabolism, oxidative phosphorylation, and others (Figures 4C,D). For the CC category, genes containing DMMSs were mainly enriched in the intracellular organelles, cytoplasm, and nucleus intracellular membrane-bound organelles. For the MF category, upregulation of m6A was notably enriched in genes involved in enzyme binding, actin binding, DNA binding, and RNA compound binding, while loss of m6A methylation was enriched in genes involved in catalytic activity, ion binding, coenzyme binding, and cofactor binding (Figures 4A,B). These results suggest that m6A has complicated roles in CI-AKI, with primary roles in metabolism, various pathways related to cell death, and oxidation.
In TvI (Figure 5), for the BP category, genes with increased methylation of m6A sites were significantly (p < 0.05) enriched in some development-associated processes, such as tissue, skeletal system, epithelium, and organ development (Figure 5A), while those with decreased methylation were highly enriched in acid metabolism processes, such as the oxoacid, organic acid, and carboxylic acid metabolic processes ( Figure 5B). Regarding pathways, the former were significantly involved in gap junctions, protein digestion and absorption, and the Wnt pathway, while the latter were involved in metabolic processes, such as the PPAR signaling pathway (Figures 5C,D). For the CC category, genes containing DMMSs induced by cisplatin-induced AKI were mainly enriched in plasma membrane, cell periphery extracellular vesicular exosome, and extracellular organelles. For the MF category, increased methylation of m6A was notably enriched in calcium ion binding, protein binding, and channel activity, while decreased methylation of m6A was enriched in enzyme binding, cell adhesion molecule binding, and identical protein binding (Figures 5A,B). These results suggest that m6A methylation is involved in berberine alleviating CI-AKI in mouse kidneys.
Finally, GO enrichment and KEGG pathway analyses were conducted on genes containing m6A sites showing opposite methylation trends in IvC and TvI (Figure 6). For the BP category, genes with DMMSs were significantly (p < 0.05) enriched in hormone secretion processes, such as regulation of hormone secretion and peptide hormone secretion ( Figure 6A). Pathway analyses demonstrated that genes with DMMSs were involved in platelet activation, complement and coagulation cascades, and gastric acid secretion (Figure 6B). For the CC category, genes containing DMMSs were mainly enriched in the cell periphery, plasma membrane, and extracellular matrix. For the MF category, genes were enriched in binding, cell adhesion molecule binding, and calcium ion binding ( Figure 6A). These results suggest that berberine might resist the nephrotoxicity of cisplatin through different pathways.
To further analyze the role of m6A in alleviating the action of berberine, genes with opposite methylation and expression trends in IvC and TvI were selected (Table 3). Overall, 9 were hypermethylated and upregulated in IvC but hypomethylated and downregulated in TvI, and 12 were hypomethylated and downregulated in IvC but hypermethylated and upregulated in TvI. GO and pathway analyses were performed to uncover the biological processes associated with these genes (Figure 7C). These genes were found to be highly enriched in complement and coagulation cascades, platelet activation, and cytokine-cytokine receptor interaction (Figure 7D). For instance, FGA/FGB/FGG genes encoding fibrinogen (Martini et al., 2019;Vilar et al., 2020) were hypermethylated and upregulated in IvC but were hypomethylated and downregulated in TvI ( Figure 8A). Slc12a1, the key gene that mediates the electroneutral movement of Na(+) and K(+) across cell membranes (Schiessl et al., 2013;Xue et al., 2014;Hao et al., 2018;Kawaguchi et al., 2018), was hypomethylated and downregulated in IvC but hypermethylated and upregulated in TvI ( Figure 8B). Havcr1, also known as Kim-1, is a well-known biomarker for kidney injury (Lippi et al., 2018;Seibert et al., 2018;Song et al., 2019;Zdziechowska et al., 2020). Cisplatin induced hypermethylation and upregulation of Havcr1, while berberine reversed these effects ( Figure 8C).

Western Blotting Analyses of Protein Expression in FGA, Slc12a1, and Havcr1
Western blotting was conducted to detect the protein levels of FGA, Slc12a1, and Havcr1 in each group. Cisplatin induced an increase in Slc12a1 protein levels and a decrease in FGA and Havcr1 protein levels. However, berberine pretreatment reversed these effects (Figure 8D).

DISCUSSION
m6A methylation is considered a reversible dynamic modification in many species. Previous studies have shown that m6A modification can regulate cellular responses to stimuli by affecting mRNA transcription, splicing, localization, translation, stability, and posttranscriptional regulation of gene expression at the RNA level (Zhou et al., 2015;Fry et al., 2017). Evidence (Anders et al., 2018;Wang Y. et al., 2019;Zhou P. et al., 2019;Li et al., 2020;Liu et al., 2020;Xu et al., 2020) also suggests a strong relationship between m6A modification and kidney disease, thereby revealing an important biological role for m6A in the regulation of kidney injury. In this study, a reliable cisplatin-induced AKI model was established in mice, and the model was tested by analyzing Scr, BUN, and kidney section images. The kidney injury resulted in dynamic m6A modifications and gene expression in the kidneys. Differentially methylated mRNAs were found to be involved in many biological pathways. m6A methylation was enriched in the CDS in all groups. In IvC, GO, and KEGG analyses of coding genes harboring DMMSs demonstrated that genes with increased methylation were primarily enriched in the pathways related to metabolic processes and cell death process, such as macromolecule metabolic processes, cellular metabolic processes, the TNF signaling pathway, and the MAPK signaling pathway, while decreases in methylation were mainly enriched in metabolic processes, oxidation, and transport, indicating that cisplatin induced complicated variation in m6A methylation. Several pathways have been thoroughly studied as factors affecting CI-AKI. However, we found that variation in metabolic processes constituted the foundation of the biological and pathological changes in CI-AKI. Different m6A methylation sites were also detected in TvI. GO, and KEGG analyses of coding genes with DMMSs revealed that the genes with increases in methylation were primarily enriched in pathways associated with tissue, the skeletal system, and epithelium development, such as gap junctions and the Wnt signaling pathway. Genes with decreases in methylation were mainly enriched in the oxoacid metabolic process and organic metabolic process. These results provide evidence of a link between m6A and berberine-mediated regulation of kidney injury, showing that berberine attenuated CI-AKI by increasing the methylation of genes associated with tissue, the skeletal system, and epithelium development. These results should prove useful for directing future research into the underlying medical value of berberine in CI-AKI.
Furthermore, 94 DMMSs showed opposite methylation trends in IvC and TvI. The genes with those DMMSs were primarily enriched in pathways associated with hormone secretion and complement and coagulation cascades, which implies that berberine probably relieves CI-AKI by directly reversing some changes in m6A methylation.
Finally, conjoint analyses of m6A modifications and gene regulation provided a broad picture of CI-AKI and the impacts of berberine pretreatment. In total, 21 genes were found to show opposite m6A methylation and expression trends in IvC and TvI. Among these candidate genes, several drew our attention. For instance, FGA/FGB/FGG are the core genes encoding fibrinogen. Cisplatin may cause endothelial injury and inflammation, which can activate coagulation cascades (Ozkok and Edelstein, 2014). Thus, FGA/FGB/FGG were hypermethylated and upregulated in CI-AKI. However, berberine pretreatment significantly alleviated these trends. This result is consistent with one previous study (Riccioni et al., 2018). Slc12a1, also known as NKCC2, is the molecular target of loop diuretics, as it is expressed on the apical membrane of the thick ascending limb of Henle epithelial cells (Hao et al., 2018). It mediates the electroneutral movement of Na (+) and K (+) , which is tightly coupled to the movement of Cl (−) across cell membranes (Schiessl et al., 2013;Xue et al., 2014;Kawaguchi et al., 2018). In our study, Slc12a1 was hypomethylated and downregulated by cisplatin. Berberine pretreatment maintained transmembrane ion exchange by hypermethylating and upregulating Slc12a1. Havcr1, also known as KIM-1, is a promising biomarker for predicting kidney injury (Lippi et al., 2018;Seibert et al., 2018;Song et al., 2019;Zdziechowska et al., 2020). It may reduce acute injury by mediating phagocytosis . In our study, it was hypermethylated and upregulated in CI-AKI. However, berberine reversed these trends. Consistently, the protein expression levels of FGA, Slc12a1, and KIM-1 were similar to the m6A methylation and gene expression levels. These results suggest that berberine has a protective effect in CI-AKI through different pathways.

CONCLUSION
We characterized the differential m6A methylome in normal, CI-AKI, and berberine-pretreated CI-AKI in mice. Our results suggest a strong association between m6A methylation and the regulation of CI-AKI. In addition, the candidate genes identified reveal the possible pathways by which berberine alleviates kidney injury. Altogether, our results provide a fundamental contribution to research into the mechanisms of and novel therapies for CI-AKI.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in (US National Center for Biotechnology Information Gene Expression Omnibus). The accession ID is GSE157261 (https://www.ncbi.nlm.nih.gov/geo/query/ acc. cgi?acc = GSE1 57261).

ETHICS STATEMENT
The animal study was reviewed and approved by The Institutional Animal Care and Use Committee of Shanghai Jiaotong University affiliated Renji Hospital.

AUTHOR CONTRIBUTIONS
ZN and XC contributed to design the research. JS, WW, and XS conducted the experiments. JS and WW analyzed and interpreted the data, and drafted the manuscript. JW and SL helped to polish the manuscript. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
We are grateful to BioDataStudio (Shanghai) for the assistance of drawing figures. We are also grateful to Drs. Chaojun Qi and Wenyan Zhou for the assistance in pathobiology.