Analysis of potential biomarkers for diabetic kidney disease based on single-cell RNA-sequencing integrated with a single-cell sequencing assay for transposase-accessible chromatin

Single-cell transcriptome analysis of the kidneys from db/db and db/m mice

We performed scATAC-seq on kidneys of 16-week-old db/db and db/m mice to characterize the accessible chromatin landscape of kidneys stimulated by hyperglycemia. Transmission electron microscopy (TEM) revealed obvious structural changes in db/db mouse kidneys compared to those in db/m mouse kidneys, including increased accumulation of the extracellular matrix, thickened basement membrane, and fusion of foot processes. (Figure 1A). A total of 42,833 cells from six kidneys obtained from three db/db and three db/m mice were processed for scRNA-seq on the 10× Genomics Chromium platform, with an average of 57,099 reads per cell. After quality control, 32,058 cells were included in the analysis (Supplementary Table 1), comprising approximately 13,906 cells from db/db mice and 18,152 cells from db/m mice, with an average of 2,096 genes per cell.

Figure 1. Single-cell transcriptome analysis in db/db mice. (A) Transmission electron micrographs of glomerulus from db/db and db/m mouse kidneys. The red arrow points towards the foot process fusion on the glomerular basement membrane. The green bar in the db/db image on the right side indicates the thickening of Glomerular basement membrane. (B) tSNE visualization of kidney cells from db/db and db/m mice, colored by cell type. DCT, distal convoluted tubule; PTS3, third segment of the proximal tubule; PTS2, second segment of the proximal tubule; PTB, proximal tubule brush; CD, collecting duct; LOH, loop of Henle; PTS1, first segment of the proximal tubule; B cell, B lymphocytes; M1, classical macrophage; T cell, T lymphocyte; M2, alternatively activated macrophage; Endo, endothelial cells; Neut, neutrophil progenitor; Podo, podocyte; NK, natural killer cell. (C) Bar plot showing the percentages of transcriptionally defined cell populations in db/db and db/m mice. The colors correspond to the cell types identified. (D) Heatmap of the mean expression of two manually selected marker genes in each cell type. Gene expression was standardized between –3 and 3 and is indicated by color intensity.

A graph-based clustering method based on the standardized expression of marker genes (Supplementary Table 2) was used to annotate the following 15 distinct populations, which were visualized in the t-distributed stochastic neighbor embedding (tSNE) space (Figure 1B): the distal convoluted tubule (DCT; expressing Slc12a3, Pvalb), third segment of the proximal tubule (PTS3; expressing Atp11a, Malat1), second segment of the proximal tubule (PTS2; expressing Cndp2, Nudt19), proximal tubule brush (PTB; expressing Hbb.bs, Hba.a2), collecting duct (CD; expressing Atp6v1g3, Atp6v1b1), loop of Henle (LOH; expressing Umod, Slc12a1), first segment of the proximal tubule (PTS1; expressing Slc5a2, Spp2), B lymphocytes (B cells; expressing Cd79a, Cd19), T lymphocytes (T cells; expressing Cd3e, Cd3g), classical macrophages (M1; expressing Lyz2, Cfh), alternatively activated macrophages (M2; expressing C1qb, Cd74), endothelia (Endo; expressing Bst2, Ccl4), neutrophil progenitors (Neut; expressing S100a9, S100a8), podocytes (Podo; expressing Cryab, S100a6), and natural killer cells (NK; expressing Gzma, Ccl5). The transcriptional signatures of these 15 cell clusters were identified. The cell numbers in each cluster in db/db and db/m mice are listed in Supplementary Table 3, and their proportions are shown in Figure 1C. Proximal tubule cells were the most abundant in all clusters, accounting for approximately 85% of total cells, whereas podocytes were the least abundant, accounting for approximately 0.3% of cells, followed by collecting duct cells. Compared with db/m mice, db/db mice had significantly decreased cell numbers in the Podo, PTS3, PTS2, and Endo clusters (P<0.01) but significantly increased cell numbers in the LOH, CD, PTB, PTS1, and DCT clusters (P<0.01, Figure 1B and Supplementary Table 3). The numbers of M1, M2, T, B, and Neut cells significantly decreased in db/db mice (P<0.01).

The top two marker genes from each cell type ranked by significance after standardization are shown in Figure 1D. We analyzed the expression patterns of specific transporters in the renal tubules. Genes encoding aquaporin 11 (Aqp11) and transporters for glucose (Slc2a2), organic cations (Slc22a1), organic anions (Slc22a8), amino acids (Slc3a1, Slc7a9, Slc1a1, and Slc6a18), uric acid (Slc22a12), sulfate (Slc13a1), phosphate (Slc34a1), and thyroid hormone (Slc16a2) were specifically expressed in the proximal tubules (Supplementary Figure 1A–1C). The expression of gene encoding aquaporin 6 (Aqp6) was specific to the CD, the expression of genes encoding sodium and chloride iron transporter (Slc12a3) and serine/threonine-protein kinase (Trpm6) were specific to the DCT, the expression of gene encoding sodium, potassium and chloride ion transporter (Slc12a1) was specific to the LOH (Supplementary Figure 1D).

Comparative analysis of renal tubule cells from db/db and db/m mice

PTS3 had the most differentially expressed genes (DEGs) between db/db and db/m mice (n = 6,284), followed by PTS1 and PTS2, whereas CD had the least DEGs (n = 140), followed by Neut (Supplementary Table 4). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that the significantly upregulated genes in the PTS3 cluster of db/db mice were mainly involved in the PPAR signaling pathway and fatty acid absorption and degradation (Figure 2A). In the tubules, such as LOH, CD, and DCT, the significantly upregulated genes in db/db mice were mainly involved in the MAPK and TGF-β signaling pathways, whereas those in Podo and Endo cells were mainly involved in the Wnt and MAPK signaling pathways (Supplementary Figure 2).

Figure 2. Comparative analysis of differentially expressed genes between db/db and db/m mice. (A) Functional enrichment analysis of the upregulated genes in the PTS3 cluster based on KEGG pathway analysis. The color corresponds to the significance values. (B) Heatmap showing the numbers of cell–cell interactions in the scRNA-seq datasets of db/m mice and db/db mice. The strength of cell–cell communication is indicated by color intensity. (C) Differential expression of Txnip and Wfdc2 in db/db and db/m mice. (D) Wfdc2 expression in the morning urine and kidneys of patients with diabetic kidney disease (DKD) and non-diabetic kidney disease (NDKD) controls, as analyzed by western blotting and immunohistochemistry. (E) Violin plots showing the expression of Klf10 and Klf11 in the kidneys of db/db and db/m mice.

Analysis of cell communication showed a significant enhancement in the kidneys of db/db mice, especially between podocytes and tubular cells (Figure 2B). Similarly, the analysis of intercellular receptor–ligand pairs revealed that the expression of ligand receptors for Sema4D-PlexnD1 and Sema4D-PlexnB1 was higher and that of ligand receptors for VEGFB_ADRB2 and VEGFB_NRP1 was significantly lower in db/db mice than in db/m mice (P<0.01, Supplementary Figure 3A).

Moreover, we found that the expression of some kidney damage markers increased in db/db mice (Figure 2C, 2E and Supplementary Figure 3B). The expression of Klf10 and Klf11 showed the most significant increase in the proximal tubules and the PTS3 cluster, respectively. On the other hand, Arg2 expression significantly increased in the LOH and CD clusters (Supplementary Figure 3B). The expression levels of Txnip and Wfdc2 were significantly increased in the proximal tubules of db/db mice compared to those of db/m mice (P<0.01, Figure 2C). As Wfdc2 is a secretory protein associated with renal fibrosis, we chose it for further analysis in DKD patients. We compared Wfdc2 expression in the morning urine of patients with DKD and NDKD controls. Western blotting showed that Wfdc2 expression was higher in the morning urine of patients with DKD than in that of controls (Figure 2D). Immunohistochemistry revealed high Wfdc2 expression in the lumen of the renal tubules and moderate Wfdc2 expression in the cells of the renal tubules, and the expression levels were higher in DKD patients than in NDKD controls (Figure 2D).

Metabolic reprograming in db/db mice

The expression levels of glucose reabsorption transporters (Slc5a1, Slc5a2, Slc2a5, Slc5a9, and Slc5a10) significantly increased in the proximal tubules of db/db mice compared with those in the proximal tubules of db/m mice (P<0.01) (Figure 3A, 3B). Reabsorption of glucose in the renal tubules is mainly completed through SGLT2 (Slc5a2) (90%) [12]. We performed immunohistochemical analysis to localize the expression of SGLT2 (Slc5a2) in the kidneys of db/db and db/m mice. SGLT2 was found to be expressed in the luminal membrane of proximal renal tubules and the urinary pole of Bowman’s capsule in db/db mice at a significantly higher level than in db/m mice (Several are the whole Bowman’s capsule) (Figure 3B).

Figure 3. Metabolic changes in the renal tubule cells of db/db and db/m mice. (A) Violin plots showing the expression of glucose transporter in db/db and db/m mice. (B) Violin plots and immunohistochemical staining results of the expression of the Slc5a2 (SGLT2) transporter in db/db and db/m mice. Scale bar, 20 μm. (C) Violin plots showing the expression of Apoe, Slca27a1, Cox6b2, and Hmgcs2 transporters in db/db and db/m mice.

Regarding lipid transport, Apoe expression level significantly increased in the proximal tubules and that of Slc27a1 significantly increased in the PTS3 cluster of db/db mice than that of db/m mice (P<0.01) (Figure 3C). Additionally, the expression of Cox6b2, involved in the oxidative phosphorylation of glucose, significantly increased in db/db mice compared with that in db/m mice (P<0.01), especially in the proximal tubules, CD, and LOH (Figure 3C). Hmgcs2, which is involved in ketogenesis, was also highly expressed in the proximal tubules, endothelial cells, and podocytes of the kidneys of db/db mice compared with that in db/m mice (P<0.01, Figure 3C). Conversely, the expression levels of the amino acid transporters Slc1a1, Slc3a1, and Slc13a1 were lower in the proximal tubules of db/db mice than in those of db/m mice (P<0.01, Supplementary Figure 3C). Transporters related to sugar and ketone metabolism were highly expressed in db/db mice; however, the expression of some key amino acid transporters decreased. These data suggest that db/db mice undergo metabolic remodeling, manifested as metabolic changes in glycolipid and amino acid levels.

Single-cell chromatin accessibility analysis in db/db and db/m mice

We sequenced 16,660 cells from the kidneys of 3 db/db and 3 db/m mice, with an average of 7,069 fragments per cell; 14,968 cells were screened for further analysis after quality control. The cells from scATAC-seq were annotated as 15 clusters under the reference of defined cell types from the scRNA-seq data, which were visualized using tSNE (Figure 4A). We identified 85,212 differentially accessible open chromatin peaks (DAPs) across the 15 cell types (Figure 4B and Supplementary Table 5). The expression of 884 cell type specific TFs which are likely to bind to regions of high chromatin accessibility across the 15 clusters was observed (Figure 4C and Supplementary Table 6).

Figure 4. Single-cell chromatin accessibility analysis of db/db and db/m mice. (A) tSNE of scATAC-seq data from db/db and db/m mice kidney cells, colored by cell type. B cell, B lymphocytes; CD, collecting duct; DCT, distal convoluted tubule; Endo, endothelial cells; LOH, loop of Henle; M1, classical macrophage; M2, alternatively activated macrophage; NK, natural killer cell; Neut, neutrophil progenitor; Podo, podocyte; PTB, proximal tubule brush; PTS1, first segment of the proximal tubule; PTS2, second segment of the proximal tubule; PTS3, third segment of the proximal tubule; T cell, T lymphocytes. (B) Cell type-specific peaks in scATAC-seq data. (C) Cell type-specific transcription factor (TF) expression based on integrative analysis of scATAC-seq and scRNA-seq data. (D) Browser tracks of marker genes from scRNA-seq data in each cell type based on scATAC-seq data. (E) Heatmap showing the expression of cell type-specific positive TFs in each cluster, with the expression level indicated by color intensity.

Browser tracks of marker genes from scRNA-seq data in each cell type were shown based on scATAC-seq data, including 6v1g3 for CD, Slc12a3 for DCT, Ccl4 for Endo, Slc12a1 for LOH, Cryab for Podo, and Slc5a2 and Atp11a for the proximal tubules (Figure 4D). We identified the cell type specific positive TFs whose expression levels were positively correlated to changes in the accessibility of their corresponding motif by integrating the scRNA-seq and scATAC-seq data (Figure 4E and Supplementary Table 7). The expression of cell type specific positive TFs in each cluster is displayed as a heatmap (Figure 4E). Proximal tubule specific positive TFs included transcription factor Maf (Maf), hepatocyte nuclear factor 4-alpha (Hnf4a), homeobox protein cut-like 1 (Cux1), bile acid receptor (Bar/Nr1h4), and hepatocyte nuclear factor 1-beta (Hnf1β). The endothelial-specific positive TFs included transcription factor 4 (Tcf4) and X-box-binding protein 1 (Xbp1), and the podocyte-specific positive TFs included forkhead box protein Q1 (Foxq1) and GATA-3 (Gata3).

Footprinting showed the predicted precise binding location of Hnf4a and Hnf1β, and tSNE was used to visualize their expression in the proximal tubules (Figure 5A–5D and Supplementary Figure 4). G6pc is a target gene of Hnf4a [13], and its expression was consistent with that of Hnf4a, mainly being expressed in the proximal tubules (Figure 5E). Rbm47 is a target gene of Hnf1β [14], and its expression was consistent with that of Hnf1β (Supplementary Figure 5E). Footprinting also showed the predicted precise binding location of Maf, Cux1, and Nr1h4, and the chromatin accessibility of the regions that Maf, Cux1, and Nr1h4 bound were visualized using tSNE (Supplementary Figure 4F–4H).

Figure 5. Specific transcription factors in the proximal tubules. (A) Footprinting of Hnf4a. (B) The chromatin accessibility of the region that Hnf4a bound was visualized using tSNE, with the accessibility level indicated by color intensity. (C) Differential motif enrichment analysis of Hnf4a: density map showing the distribution of motif deviation coefficient values. (D, E) Expression of Hnf4a and G6pc (target gene of Hnf4a), with the expression level indicated by color intensity. (F) Peak distribution, footprinting, and expression of Hlf in the PTS3 cluster of db/db and db/m mice. (G) Peak distribution, footprinting, and expression of Dbp in the PTS3 cluster of db/db and db/m mice. (H) Peak distribution, footprinting, and expression of Tef in the PTS3 cluster of db/db and db/m mice.

Circadian rhythm-associated TFs in db/db renal tubular cells

Hepatic leukemia factor (Hlf), D site-binding protein (Dbp), and thyrotroph embryonic factor (Tef) were found to bind to the regions with high chromatin accessibility in the PTS3 cluster of db/db mice more than that in db/m mice (P<0.01, Figure 5F–5H). Hlf, Dbp, and Tef belong to the PAR-domain basic leucine zipper (PAR bZip) family and serve as the output executors of the biological clock [15, 16].

Hypoxia-inducible factor 1-alpha (Hif1α) and Kruppel-like factor 9 (Klf9) showed higher chromatin accessibility of their promoter regions in the PTS3 cluster of db/db mice than in those of db/m mice (P<0.01, Figure 6A, 6B). Nuclear respiratory factor 1 (Nrf1) and pre-B-cell leukemia transcription factor 1 (Pbx1) showed lower chromatin accessibility of their promoter regions in the PTS1 cluster of db/db mice than in those of db/m mice (P<0.01). Pbx1 is necessary for renal development and is involved in the activation of the nephrin promoter in the proximal tubules [17], while Nrf1 acts as a cytoprotective factor in cells and plays a regulatory role in diabetes by modulating the level of glutathione [18] (Figure 6C, 6D).

Figure 6. Kidney injury-related transcription factor (TF) analysis in the PTS1 and PTS3 clusters. (A) Peak distribution, footprinting, and expression of Hif1a in the PTS3 cluster of db/db and db/m mice. (B) Peak distribution, footprinting, and expression of Klf9 in the PTS3 cluster of db/db and db/m mice. (C) Peak distribution, footprinting, and expression of Pbx1 in the PTS1 cluster of db/db and db/m mice. (D) Peak distribution, footprinting, and expression of Nrf1 in the PTS1 cluster of db/db and db/m mice.

Model animals

All animal experiments were performed in accordance with the experimental guidelines of the Experimental Animal Ethics Committee of Zhengzhou University (approval number ZZU-LAC20201013). Three 16-week-old male db/db mice and three 16-week-old male db/m mice were purchased from the Jiangsu Jicui Yaokang Biotechnology Company. The baseline information of db/db and db/m mice is listed in Supplementary Table 8.

Specimens from DKD patients and NDKD controls

The DKD urine and kidney specimens were obtained from hospitalized patients who were diagnosed with Diabetic Nephropathy based on the pathology of renal biopsies in the First Affiliated Hospital of Zhengzhou University (n=3). The NDKD kidney specimens were collected from patients undergoing surgical treatment for renal carcinoma (n=3). Specimens used in this study were collected 5cm from the tumor. This study was approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University.

Preparation of single-cell suspensions

Mouse kidney tissues were placed in a sterile Petri dish with 10 mL precooled (4° C) 1× Dulbecco’s phosphate-buffered saline (PBS, Thermo Fisher Scientific, Waltham, MA, USA). After washing twice, the tissues were cut to obtain 1- to 2-mm³ pieces, placed in fresh dishes with lysis buffer that contains 0.25% trypsin (Thermo Fisher Scientific) and 10 μg/mL DNase I (Sigma-Aldrich, St. Louis, MO, USA), and lysed in a 50-mL centrifuge tube in a water bath at 37° C. An equal volume of the complete medium (Thermo Fisher Scientific) was added to stop the lysis. The buffer was passed through a pre-wetted 70-μm filter (Millipore, Burlington, MA, USA), and the filtrate was collected in a fresh 50-mL centrifuge tube. The filter was rinsed with 10 mL complete medium, whereas the filtrate was collected in the same centrifuge tube and centrifuged at 100 ×g for 5 min at 25° C. The supernatant was discarded, and 5 mL 1× Dulbecco’s PBS containing 2% fetal bovine serum (FBS, Gibco, Grand Island, NY, USA) was added to resuspend the cells. The suspension was centrifuged at 100 ×g for 5 min, the supernatant was discarded, and 3 mL precooled (4° C) 1× erythrocyte lysis buffer (Thermo Fisher Scientific) was added, followed by gentle blowing and mixing. The suspension was incubated at 25° C for 3 min and centrifuged at 300 ×g for 5 min, and the supernatant was discarded. Next, 5 mL precooled 1× Dulbecco’s PBS containing 2% FBS was added, followed by gentle blowing and mixing; this mixture was centrifuged at 300 ×g for 5 min. These steps were repeated until a single-cell suspension was obtained. The suspension was placed on ice for cell counting.

scRNA-seq

scRNA-seq was performed using a 10× Genomics platform. The single-cell suspension, 10× barcode gel magnetic beads, and oil were added to the different chambers of chromium chip A (Chromium™ Chip A Single Cell Kit, 120236, 10× Genomics, Pleasanton, CA, USA) to generate single-cell gel beads-in-emulsion (GEM) through the microfluidic “double cross” system. The magnetic beads on the gel contained 30-nt oligo primers to reverse-transcribe poly-A RNA into cDNA chains with barcode and unique molecular identifier information. The purified cDNA chains were amplified via PCR and quantified using Qubit. Fragment size was detected using an Agilent 2100 system (Agilent Technologies, Santa Clara, CA, USA).

Optimal fragments obtained by enzyme digestion were screened using magnetic beads (SPRIselect Reagent Kit, B23318, Beckman Coulter, Pasadena, CA, USA). A cDNA library containing the P5 and P7 linkers was constructed using PCR after terminal repair and the addition of A and adapter-ligated Read2 sequencing primers (Chromium i7 Multiplex Kit, 120262, 10× Genomics; Chromium™ Single Cell 5′ GEM, Library and Gel Bead Kit v3, 1000006, 10× Genomics).

The library was purified using magnetic beads, the concentration of the fragments in the library was detected using Qubit, and the fragment size was measured using Agilent 2100. T and B cell receptor libraries were constructed similarly (Chromium™ Single Cell V(D)J Enrichment kit, 1000005 and 000016, 10× Genomics). The library was sequenced and analyzed using the Illumina sequencing platform (Illumina, San Diego, CA, USA). Sequencing was controlled using the data collection software provided by Illumina, and real-time data analysis was performed.

scATAC-seq

scATAC-seq was performed using the 10× Chromium™ Single Cell ATAC platform. The cell suspension was added to a 1.5-mL centrifuge tube, followed by centrifugation at 300 ×g (4° C) for 5 min. The supernatant was discarded, and lysis buffer was added to the pellet, followed by gentle blowing to resuspend the nuclei and incubation on ice for 15 min. The proportion of live cells was monitored using a fluorescent cell counter until it dropped to 0%–5%. Finally, the nuclei were washed and counted using a cell counter.

Nuclei were harvested and immediately analyzed on the 10× Chromium™ Single Cell ATAC platform. The nuclear suspension, 10× barcode gel magnetic beads, and oil were added to the different chambers of chromium chip H (Chromium™ Chip H Single Cell Kit, 1000161, 10× Genomics) to generate single-cell GEM through the microfluidic system. The transposase with the DNA adapter came into contact with the DNA where chromatin was loosely wound, resulting in a transposition reaction at these sites. As a result, a DNA fragment was cut from the chromatin and formed a free DNA fragment, which was then bound by the DNA adapter. The Read 1N sequence attached to the DNA label strand on the gel bead was complementary to the linker sequence on the DNA fragment and was annealed and adhered. The GEM was then transferred into the thermal cycler. Under the action of polymerase, the DNA fragments were extended to 10x the Barcode sequences, connected to the P5 sequences, used for library construction (Chromium™ Single Cell ATAC GEM, Library and Gel Bead Kit v1.1), and sequenced.

Qubit was used to measure the concentration, and the Agilent 4200 system was used to detect the fragment size. The library was sequenced and analyzed using the Illumina sequencing platform. Sequencing was controlled using the data collection software provided by Illumina, and real-time data analysis was conducted.

Data processing and quality control

Cellranger (cell range: 3.0.1) was used to map Raw FASTQ files to the mm10 reference genome (Ensembl GRCm38,93), and the scRNA-seq data were converted from the fastq type to a cell expression matrix, retaining cells with >200 expressed genes. The effects on the cell cycle were negligible and could not be detected in this study.

Dimension reduction, batch effect correction, and clustering

The top 2000 genes with the largest variation were identified based on the average value and dispersion (variance/average) of all genes in each sample. We used the combination of canonical correlation analysis and mutual nearest neighbors to correct for batch effects. Principal component analysis was performed for dimension reduction, and the first 20 dimensions were used for subsequent analysis. Normalized data were clustered and analyzed using the Louvain algorithm (resolution: 0.5) and visualized using tSNE.

Cell-type annotation and marker genes

The cell type was annotated manually using singleR and marker gene software referring to the marker gene list (Supplementary Table 2). The marker genes of all clusters were analyzed using the Wilcoxon algorithm (P < 0.05). The genes in each cluster that showed high specific expression (log-fold change [FC] > 0.25) and were expressed in at least 20% of cells were selected as significant marker genes of the corresponding cluster, displayed using R pheatmap, and visualized using tSNE.

DEG and functional enrichment analyses

The DEGs between different clusters and groups were analyzed using the FindMarkers function in Seurat, followed by KEGG functional enrichment analysis. Gene expression profiles were visualized using ggplot, and heatmaps were constructed using the R pheatmap package.

Quality control, doublet filtration, and dimension reduction

We used CellRanger to convert the scATAC-seq data from the six mouse samples from the FASTQ format to the fragment format. We calculated the transcription start site (TSS) enrichment score of each cell (average accessibility in the 50-bp region centered on each single-base TSS position divided by the average accessibility [±1900–2000 bp] of the TSS flank position). Cells with a TSS enrichment score of >4 and >3000 fragments size were eliminated (Supplementary Figure 5). ArchR was used to ascertain whether the cell is a “double” (i.e., two different genotypes within a cell), and double scores were allotted to each cell accordingly. Double cells were filtered at a cut-off rate of 5%. ArchR was used to reduce the hierarchical dimensionality, and harmony was used to correct the sample batch effects. tSNE was used for data visualization.

Clustering and annotation

Gene scores were evaluated based on the open chromatin data obtained after scATAC-seq and calculated using the distance-weighted openness model. The gene score represents the extent to which gene expression is predicted by the regulatory elements near the gene. First, the data were integrated by mapping the gene score matrices of the scATAC-seq data to those of the scRNA-seq data [41]. Second, the Seurat clustering method was used to annotate the cell types in the scATAC-seq data based on the defined clusters of the scRNA-seq data. ArchR was used to identify marker genes according to the gene scores. Significant marker genes of each cluster were identified (false-discovery rate ≤ 0.01 and log2 FC ≥ 1.25). The gene scores of each cell were embedded into tSNE to visualize the marker genes from the scATAC-seq data.

Peak calling, marker peak identification, and motif enrichment

ArchR was used for peak calling of the scATAC data, which mainly involved the construction of pseudo bulk replicates (i.e., data that imitate bulk ATAC-seq after all single cells were merged) from multiple pseudo bulk samples in each cluster. Marker peak analysis of each cluster was performed similar to cluster marker gene identification. After obtaining the marker peaks specifically enriched in each cluster [42], the getMarkerFeatures function in ArchR was used for motif enrichment analysis to predict the TFs associated with transcription-binding events and open chromatin sites [43]. The browser tracks of these TFs were drawn.

Footprint analysis, distribution of the chromVAR deviation score, and positive TF regulators

TF footprints can predict the precise binding position of TFs at specific locations. We used ArchR to analyze the Tn5 deviation footprints. Due to the deviation of Tn5 transposase insertion sequences, the distribution of the chromVAR deviation scores of TFs was assessed using chromVAR. Normalization was performed to remove the bias associated with Tn5 from the coverage signal. ArchR was used to identify the positive TF regulators in the scATAC data. We correlated the chromVAR deviation Z scores of the TF motifs with the expression of TF genes from low-overlap aggregates and selected the TFs whose motif accessibility was positively correlated with TF expression. Target genes specific to the TFs of interest were predicted using the TRRUST v2 and Cistrome Data Browser online database [44, 45].

Immunohistochemical analysis

First, the paraffin sections are dewaxed and hydrated by heating at 60° C for 30 min, treatment with xylene for 15 min, and treatment with a gradual ethanol series of 100%, 95%, 85% and 75% ethanol for 5 min each. The sections then underwent antigen retrieval using the microwave oven repair method. Sections were incubated with 3% H₂O₂ to inactivate endogenous peroxidases and then blocked with 2% goat serum for 30 min. Subsequently, the sections were incubated overnight at 4° C with primary antibodies (SGLT2 antibody, BiCell Scientific, BC-20802; Wfdc antibody, Sigma-Aldrich, HPA042302). Color rendering was performed using diaminobenzidine (DAB), and the sections were photographed.

Western blotting

The morning urine of DKD patients and healthy controls were collected and filtered using a 0.2-μm filter and then concentrated 10-fold using the Centrifugal Filter Unit (Merck Millipore, UFC9010). The samples were run on a 15% SDS-PAGE gel at 120 V for 40 min. The gel was transferred to a methanol-activated PVDF membrane and run at 250 mA for 30 min and then blocked with 5% skimmed milk. The PVDF membrane was then incubated at 4° C overnight with the primary antibody (WFDC2, Servicebio, Wuhan, China, GB111377). Subsequently, a secondary antibody was applied, followed by incubation, washing, and imaging.

Transmission electron microscopy (TEM)

The kidney specimens were fixed in glutaraldeyde immediately after obtained from db/db and db/m mice. Ultrathin sections were cut with an ultramicrotome, stained with 2% (wt/vol) uranyl acetate and lead citrate, and examined with a JEOL JEM-1400 Plus transmission electron microscope.

Statistical analysis

All statistical analyses were performed using R software and corrected using the false-discovery rate, P-value < 0.05 was considered statistically significant.

Data availability statement

Single cell sequencing data in FASTQ format are available at the Genome Sequence Archive (GSA) for Human at the BIG Data Center, Beijing Institute of Genomics (accession number PRJCA014746).