The db/db (BKS.Cg-Dock7m^{Lepr +/+ db/J}) and nondiabetic control (C57BL/6J) male mice at the age of 8 weeks were purchased from the Institute of Model Animal Research of Nanjing University, Nanjing, China. The mice were housed in the specific pathogen-free barrier environment of the Animal Experimentation Center, Xuanwu Campus, China Pharmaceutical University (CPU), and given adequate sterile drinking water and standard chow. The temperature was 25 °C ± 2 °C and the relative humidity was 40%–70% in the animal room. All mice were acclimatized for 1 week before the experiments were conducted. All experiments with the animals were approved by the Animal Ethics Committee of CPU and conducted according to the relevant experimental regulations.

As we previously reported (Diao et al., 2023; Yu et al., 2023b), DN in db/db mice was identified based on the presence of body glucose (BG) > 16.7 mmol/L and the urine micro-albuminuria to creatinine ratio (UACR) > 200 mg/g on two consecutive tests. The db/db mice with high urinary proteinuria were then randomly divided into HKC and DN groups. In HKC group, the suspension of HKC (0.84 g/kg/day, dissolved in distilled water) was freshly prepared and intragastrical administered to the db/db mice, while the equivalent distilled water was used simultaneously in DN group. The treatment of HKC in DN was carried out once a day for 4 weeks. HKC was purchased from Suzhong Pharmaceutical Group, Co., Ltd., Taizhou, China. The db/db mice with high BG but low UACR (db/db mice with no or minimal proteinuria) were grouped as T2D. BG (mmol/L), body weight (BW) (g), dietary intake (g), water intake (mL), and urine output (mL) were examined every week. By using a metabolic cage (Fengshi Laboratory Animal Equipment Co., Ltd., Suzhou, China), urine samples were collected for 12 h. Urine from db/db mice was collected weekly to examine UACR and urine albumin excretion rate (UAER). Quantitative ELISA kits (Elabscience, Wuhan, China) were used to measure albuminuria (MAU) and creatinine (Cr) levels in urine samples, and to calculate UACR (MAU/Cr) and UAER (MAU/12 h) values. After 4 weeks of HKC treatment, experimental mice were euthanized by decapitation following intraperitoneal administration of 30 mg/kg sodium pentobarbital. Kidney tissue samples were harvested and frozen in liquid nitrogen for subsequent experiments.

The kidneys were removed by cardiac perfusion with phosphate buffered saline and placed in formaldehyde tissue fixative. The fixed kidney tissues were embedded in paraffin and the blocks were sectioned at 4 μm with HistoCore Bio-Cutter (Leica Biosystem, Germany). The sections were then stained with Hematoxylin and Eosin (H&E) (BASO, Zhuhai, China) and/or Periodic Acid-Schiff (PAS) staining solution (Aifang Biological, Hunan, China) according to the standard procedures and finally mounted on a CX23 light microscope (Olympus, Japan) for analysis. As we have previously reported (Wu C. et al., 2022), in the present study, H&E and PAS staining sections from each animal were analyzed for semi-quantification of glomerular area, ratio of vacuolar and staining area by using the software of Image-Pro Plus (version 6.0.0260).

Genomic DNA was extracted from the kidney tissues by using the animal genomic DNA kit (E.Z.N.A.® Tissue DNA Kit, Omega Bio-Tek, United States). DNA integrity of extractions was inspected by using 1% agarose gel electrophoresis. The bisulfite transformation and DNA fragment purification were done by using a transformation kit (EZ DNA Methylation-Gold Kit, Zymo Research, United States). A high-throughput sequencing analysis was subsequently carried out with Illumina HiSeq 2500. After removing unknown nucleotides and low-quality read lengths from the original read lengths, the trimmed clean reads were detected.

For whole genome bisulfite sequencing (WGBS) analysis, the raw Fastq files obtained were trimmed by Trimmomatic (v0.39–2) under the conditions of allowing up to two base mismatches, a sliding window size of 4, and an average quality threshold of 15, and clean reads were generated by discarding unpaired reads and paired reads with a final length of less than 75 bp. The clean reads were compared to the mouse genome (mm10) by Bismark (v0.24.0) to generate bam files and extract methylation information. After removing loci with zero coverage, the following steps were performed using MethylKit (v1.22.0): methylation levels in the whole genome were calculated using the sliding window method (2 kb), as well as the Pearson correlation coefficients of all the samples (Akalin et al., 2012) to check the reproducibility within the administered groups and compare methylation differences between groups; the “methylation level” was determined by the fraction of methylated cytosines, i.e., the proportion of methylated cytosines to all cytosine sites in the region; the resulting clean reads were mapped to the mouse reference genome by using the Bismark software (v2.90); and the R package Methylkit (Wu et al., 2021) was used to estimate the methylation of CpG sites, promoter regions, CpG island region and gene annotation methylation status and ratios, as well as to identify differentially methylated regions (DMR) and differentially methylated site (DMS) between groups; 150 bp sliding-window regions with p-value <0.05 (with a step of 50 bp to satisfy an average coverage of reads greater than 10) were taken as the final DMR; base sites with q-values <0.01 and with a ratio of differences in methylation levels between base site groups of more than 25% of the base loci were DMS. DMR and DMS that overlapped with gene bodies or 2 kb regions upstream or downstream of the body region were considered as differentially methylated genes (DMGs), and genes with meth. diff >0 were categorized as hypermethylated genes (hyper genes), and genes with meth. diff <0 as hypomethylated genes (hypo genes).

RNAs were extracted from crushed kidney tissues with Trizol (Invitrogen, Carlsbad, United States), and total RNA integrity was detected by 1.2% agarose gel electrophoresis (meeting 28S rRNA/18S rRNA = 2.0). RNA concentration and purity were determined by Nanodrop 2000 UV spectrophotometer (A260/A230 = A260/A280 and >1.8), and RNA integrity was detected by Agilent 2100 (RIN value >9.0). Total RNA that met the criteria was reverse transcribed to cDNA, and ultrasonically sheared for Illumina library preparation (Illumina Truseq RNA sample prep Kit, Illumina, United States). After the library quality testing, the qualified libraries were put onto the Illumina Hiseq platform for PE150 sequencing.

For RNA transcriptome sequencing (RNA-Seq), the raw reads were trimmed as same as the procedure in the WGBS section. Trimmed data were aligned to the mouse genome (mm10) by Hisat2 (v2.2.1) to generate sam and bam files, the bam files were transformed to the gene expression count expression matrix by featurecount (v2.0.1). The expression matrices were analyzed for differential expression using DESeq2 software (v3.17), and genes with |log2FC| ≥ 0.5 and p-value <0.05 were identified as differentially expressed genes (DEGs). Among these, the genes with log2FC > 0 were categorized as upregulated genes (Up genes), while genes with log2FC < 0 were classified as downregulated genes (Down genes). Finally, the bam data were processed through rmats (v4.1.2) to output alternative 3′splice sites (A3SS), alternative 5′-splice sites (A5SS), mutually exclusive exons (MXE), retained intron (RI), and exon skipping (SE) types of variable shear events and generate sashimi-plot.

Conversion of RNA (removal of DNA) to cDNA using HiScript® III 1st Strand cDNA Synthesis Kit (Vazyme, Nanjing, China), three to four samples from each group were chosen for the experiment and RT-qPCR was carried out using Bio-Rad CFX connect™ Real-Time PCR Detection system (Bio-RAD, Singapore) and ChamQ SYBR Color qPCR Master Mix (Vazyme, Nanjing, China). The reaction system configuration and reaction conditions were performed according to the instructions. The primers were designed by the NCBI Primer-BLAST website and synthesized by Shanghai Biotech Co., Ltd., Shanghai, China (Supplementary Table S1).

To explore the function of the genes, Gene ontology (GO) and Tokyo Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using clusterProfiler (Wu et al., 2021). The pathways with P-value <0.05 were considered significantly enriched.

All quantitative results were expressed as mean ± standard error of the mean (SEM). All statistical analyses and graphics were performed using SPSS 22.0 software (SPSS, Chicago, IL) and R software (v4.3.1). Differences between groups were compared using unpaired Student’s t-test or one-way ANOVA test. The level of significance was presented as ∗ P < 0.05 and ∗∗ P < 0.01. Correlation analysis was performed by the Pearson correlation coefficient method. P < 0.05 were considered statistically significant.

We initially monitored various physiological indices in all studied animals to evaluate the reno-protective efficacy of HKC in DN. The data indicated that the UACR levels in DN group were expectedly higher compared to those in T2D group (Figure 1A). After consecutive 4 weeks of HKC administration, however, UACR and UAER levels in HKC group were significantly reduced (Figure 1A; Supplementary Figure S1A). BG and BW levels among T2D, DN, and HKC groups were analyzed, and no statistically significant change was observed (Figures 1B,C). Furthermore, H&E and PAS staining optical photographs of kidneys demonstrated that the thickened glomerular basement membrane and diffused hyperplasia, glomerular atrophy, and hyaline capillaropathy were presented in DN group while these damages were observed to be decreased in the HKC group (Figure 1D; Supplementary Figure S1B,C).

We then investigated the DNA methylation status in DN, and its changes after HKC administration by using WGBS. Based on the data of DNA methylation levels (the fraction of methylated cytosine) among T2D, DN, and HKC groups and the principal component and Pearson correlation analyses, the stabilized consistency of DNA methylation levels within each group and the obvious heterogeneity between the groups were found (Supplementary Figure S3A,B). Overall, most of the cytosines in the CpG sites were methylated. Compared to T2D, the DNA methylated sites in DN were annotated in intergenic and intron regions. The comparison of HKC with DN, however, showed that the proportion of DNA methylated sites annotated to exon and promoter regions was increased (Figure 2A), which was consistent with the percentage of methylated cytosines annotated to hypomethylated and hypermethylated genes (Figure 2B). Interestingly, there were multiple distinct methylation patterns for the multi-site in the same gene. Most of the pivotal genes in DN had three distinct types of methylation sites, i.e. CpG, CHG, and CHH while the highly epigenetic susceptibility genes after HKC administration exposure had CpG and CHG methylation patterns but not CHH (Supplementary Figure S1C).

Tet2 encode the demethylation enzymes, respectively, and play a crucial regulatory role in DNA methylation (Williams et al., 2011). In the current study, the expression of Tet2 genes was found to be increased in DN. After HKC treatment, the expression was decreased (Figure 2C). Furthermore, the aberrant DNA methylation patterns in DN group were found to be exacerbated with the progression of DN, while HKC may reverse the Tet2 expression and resulting in affecting the abnormal methylation patterns.

We further integrated WGBS and RNA-seq data to analyze gene expression levels in two comparisons: DN vs. T2D and HKC vs. DN. The different gene expression patterns of DN vs. T2D and HKC vs. DN were represented in the utilizing volcano plots (Figures 3A,B). The initiation of DN was accompanied by downregulation of massive genes (n = 6,494). The number of downregulated DEGs (n = 855) was greater than the number of upregulated DEGs (n = 785) after HKC treatment. In Figure 3C, a Circos plot showed that the DMRs in the context of CpG annotated chromosomal distribution were substantially identical, and the changes of methylation regions were evenly distributed across the genome. Overall, the development of DN and HKC administration was accompanied by altered gene methylation patterns as well as fluctuations in gene expression.

A total of 3,233 DMGs and 11,528 DEGs were identified in DN vs. T2D, along with 438 DMGs and 1,640 DEGs in HKC vs. DN. The probable biological roles of the differential genes were then elucidated through GO and KEGG pathway enrichment analyses. We found that DEGs in T2D group were mainly enriched in renal and neural development, ion transmembrane transport-relative, PI3K-Akt signaling, focal adhesion, Rap1 signaling, amino acid and fatty acid metabolism pathways (Figure 4A left). Meanwhile, DMGs were mostly involved in cell adhesion and metastasis, renal unit and glomerular formation, PI3K-Akt signaling pathway, MAPK, Ras, and Rap1 signaling pathways (Figure 4A right), in which the cross-talk between PI3K-Akt and MAPK pathway was associated with renal interstitial fibrosis (Zhang et al., 2021). HKC administration-driven DEGs were mainly enriched in biological activities such as immune regulation, viral defense, cytotoxicity, and endocytosis, as well as related to neuroactive receptor-ligand interactions and cytokine interaction pathways, and its DMGs were mainly involved in biological processes such as neural synaptic signaling, ion transmembrane transport, and Notch signaling pathway (Figure 4B). Literature records have demonstrated that Notch signaling regulates nephron number and segmentation (Bonegio and Susztak, 2012), impacting albuminuria, glomerulosclerosis, renal function, and susceptibility to renal disease (Murea et al., 2010). Inhibiting the Notch signaling system is thought to be a novel therapeutic method for DN (Lin et al., 2010).

Finally, we combined RNA-Seq and WGBS data to identify genes with both differential expression and differential methylation, revealing DN-associated methylation candidate genes that may be causal in DN and contribute to the epigenetic mechanisms underlying DN. The results of DN vs. T2D yielded 1,029 genes with one or more methylated sites and significantly altered expression with the development of DN (Figure 5A), and further investigation revealed that 577 genes were downregulated in expression due to hypermethylated sites and 6 hypomethylated upregulated genes (Figure 5B). Supplementary Table S1 shows the expression and methylation sites of 25 crucial genes involved in the methylation status in DN. Pcdh15, Mdga2, Fmo6, and Sp140 demonstrated a strong positive correlation between methylation sites and gene expression (Supplementary Figure S2). We found that hypermethylated downregulated genes were associated with kidney development and ion transport, whereas hypermethylated upregulated genes were mostly connected with biological processes including cell adhesion and synapse formation. We also discovered that numerous intersecting genes were implicated in the neuroactive ligand-receptor interaction, PI3K-Akt signaling network, cAMP signaling pathway, cell adhesion, and Ras signaling pathway. Simultaneously, there is crosstalk among various pathways (Figure 5C), forming a regulatory network strongly linked to the development of DN.

We identified Ntrk2 as a potential epigenetic susceptibility gene for DN (Supplementary Table S1). A positional candidate genetic study has demonstrated that Ntrk2 is a glomerular filtration rate variant gene and its biofunction is involved in the MAPK pathway (Thameem et al., 2015). Moreover, Rbms3 (Klumpers et al., 2022), Morc1 (Wang et al., 2021), Cd300lf (Zhang et al., 2022), and Arel1 (Rydbirk et al., 2020) are reported to be associated with epigenetic regulation in kidney diseases, while Aim2 inflammasome has potential pathogenic effects in kidney diseases, including podocyte damage and kidney inflammation (Chung et al., 2021). Podocyte injury may directly contribute to proteinuria, and the mouse podocyte autophagy regulatory protein in which Gpr137b may play an important role (Hu et al., 2020).

We further analyzed the genes with changes in methylation and transcript abundance after HKC administration, as shown in Figure 6A and Supplementary Table S2. In total, we identified 28 genes with altered expression due to DNA methylation changes. We comparatively investigated the DNA methylation and transcriptomic expression levels in db/db mice with and without HKC treatment to explore the pharmacological targets of A. manihot (L.). As shown in Figure 6B, there are 12 key genes, including Nectin1, Lars2, Zc3h7a, Cdk8, Slc16a2, Myom2, Slc22a23, Ptprd, Pde4d, Pisd-ps3, Sp140, and Sp110 screened by HKC vs. DN. We ultimately identified these 12 genes as epigenetic target genes for the treatment of DN by HKC.

DNA methylation can not only lead to changes in expression, but also cause alternative splicing (AS) during transcription (Jaenisch and Bird, 2003). We identified variable AS of each group as described above using rMATS, combining with the number of CpG methylated sites of genes, we found that HKC administration exposure not only altered the number of methylated sites but also increased the possibility and variability of AS in Zc3h7a, Pde4d, Cdk8, and Pisd-ps3, suggesting the high apparent susceptibility of these four genes might be related to the therapeutic mechanism of HKC (Supplementary Figure S4).