Comprehensive Diagnostics of Diabetic Nephropathy by Transcriptome RNA Sequencing

Introduction

Diabetes mellitus (DM) is a serious problem affecting human health at present. The International Diabetes Federation estimates that the number of adults with diabetes in the world will rise to 643 million by 2030 and there will be more than 783 million adults living with diabetes in 2045¹. More than that, the rising number of patients with diabetes has brought a huge economic burden on the medical and health care systems worldwide.²

In addition to having great harm to the endocrine system, vascular-related complications caused by diabetes can also lead to multi-system damage.³ Whereas microangiopathy in the kidney usually causes chronic kidney disease (CKD), which is often called DN⁴ and is also one of the most important and serious microvascular complications of diabetes.⁵ DN is the primary cause of end-stage renal disease (ESRD).⁶ Long-term hyperglycemia could lead to the loss of kidney function, characterized by glomerular hypertrophy, albuminuria, decreased estimated glomerular filtration rate, and renal fibrosis.⁷ It is generally believed that the pathogenesis of DN may be related to oxidative stress, activation of the RAAS system, inflammation and so on.^8,9 The early stage of DN is mainly shown in glomerular hyperfiltration, but with the progression of the disease, microalbuminuria, obvious albuminuria and progressive decrease of eGFR will be found, resulting in ESRD.¹⁰ Currently, with evaluating the longtime of diabetes or the presence of diabetic retinopathy exception, the diagnosis of DN is mainly judged by persistent albuminuria and declining eGFR.^11,12 Although these biomarkers are widely used clinically, they are inevitably affected by age, gender, heredity, and other factors. In addition, a single biomarker does not have the precise ability to predict disease.¹³ For example, some studies have shown that microalbuminuria was considered to be the prime indicator of DN though, it was not as accurate as we expected.^14–16 Given this, it is very necessary to keep seeking new biomarkers to explore the effective ways of diagnosis and treatment of DN.

In this study, three hub genes (DDX58, SAMD9L, TLR6) related to DN were obtained by bioinformatics analysis based on transcriptome sequencing data. These three genes could distinguish between DN samples and normal samples greatly, suggesting that these genes were expected to be novel DN biomarkers and potential adjuncts for clinical diagnosis.

Materials and Methods

Patient Preparation

In this study, 3 DN patients and 3 healthy volunteers in the Second Hospital Affiliated to Kunming Medical University were selected in 2020 and their blood was taken as the test samples. All the blood samples were obtained separately for the needs of the study. Previously, we ensured that each subject was fully informed the purpose and significance of the study, and the participation was entirely voluntary. All the patients were clinically diagnosed with DN and the inclusion criteria were based on the clinical diagnostic criteria with reference to the expert consensus on clinical diagnosis of DN in Chinese adults:¹⁷ the history of diabetes for more than 10 years; persistent albuminuria(>300 mg/day) or a progressive decrease in GFR; simultaneously existing diabetic retinopathy (DR). Furthermore, routine blood-examination results from the 3 healthy volunteers were within the normal range. 2mL Trizol reagent (ThermoFisher Scientific, China) was mixed with per milliliter peripheral blood and then stored at −80°C.

The study was approved by the Medical Ethics Committee of the Second Hospital Affiliated to Kunming Medical University, following the guidelines outlined in the Declaration of Helsinki. Because of the low risk to the subjects, the waiver of informed consent for the study was also agreed by the Second Hospital Affiliated to Kunming Medical University.

RNA Extraction and Qualification

Trizol (Thermo Fisher Scientific) was used to extract RNA from the samples respectively, according to the instructions. The RNA concentration and integrity were evaluated using Nanodrop 2000 spectrophotometer (ThermoFisher Scientific, USA) and Agilent Bioanalyzer 2100 (Agilent Technologies, USA).

Library Construction

Ribo-off rRNA Depletion Kit (Vazyme, China) was used to remove rRNA from each sample, and then VAHTSTM RNA Clean Beads (Vazyme, China) was used to purify rRNA-depleted RNA. According to the manufacturer’s instructions of VAHTS Universal V8 RNA-seq Library Prep Kit for Illumina (Vazyme, China), the first-strand and the second-strand of cDNA were synthesized. Then, the obtained cDNA was end-repaired, 3′adenylated. After that, an optional NEBNext adaptor was used for adaptor ligation and the product was cleaned up next. Subsequently, depending on the reaction conditions of polymerase chain reaction (PCR), the product was amplified by PCR using SYBR Green qPCR Master Mix (Servicebio, China) and purified similarly. Finally, the quality of the library was assessed by Qsep-400 (Bioptic, China) and the qualified library was sequenced on the Illumina novaseq 6000 platforms (Illumina, USA) using NovaSeq 6000 S4 Reagent Kit (Illumina, USA).

Differential Analysis

The limma package was applied to indicate DE-mRNAs and -lncRNAs between normal and DN samples, with screening criteria of P < 0.05 and |log₂ fold change (FD)| > 1.

Construction of the ceRNA Network

The ceRNA network was constructed based on DE-mRNAs and -lncRNAs. Human miRNA sequences were obtained by miRbase. In miranda software, the binding score was set to 180 to predict the relationship pairs among miRNAs, DE-lncRNAs, and DE-mRNAs. In this study, we only kept the relationship pairs with consistent expression trends, ie, up-regulated DE-mRNA-miRNA-up-regulated DE-lncRNA and down-regulated DE-mRNA-miRNA-down-regulated DE-lncRNA.

Construction of PPI Networks and Identification of Hub Genes

The Protein-Protein Interaction (PPI) network was created via the STRING website based on the mRNAs in the ceRNA network obtained above and visualized in Cytoscape. For the hub genes, the top10 mRNAs acquired with the CytoHubba plugin of Cytoscape were used to define them.

GO and KEGG Analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed in the clusterProfiler package to reveal the potential functions of hub genes. The GO system consists of three components, biological processes (BP), cellular components (CC), and molecular functions (MF). Terms with P < 0.05 were considered significantly enriched. Unfortunately, hub genes did not appear to be enriched in the CC category.

Construction and Evaluation of Diagnostic Model

We undertook LASSO regression analysis of 10 hub genes using the glmnet software package to derive the optimal variables to construct a diagnostic model for DN. ROC curves were then applied to assess the ability of the diagnostic model and individual variables to distinguish DN from normal samples.

Single-Gene GSEA

Single-gene GSEA was employed to elucidate the signaling pathways involved in diagnostic biomarkers. In this study, we calculated the correlation of each diagnostic biomarker with all other genes separately and ranked all genes from highest to lowest correlation. The ranked genes were used as the set of genes to be tested, the KEGG signaling pathway was used as the pre-defined gene set, and finally, the enrichment of the KEGG signaling pathway in the corresponding set of genes to be tested was examined.

Isolation of Peripheral Blood Mononuclear Cells (PBMC)

We collected peripheral blood samples from 10 DN patients and 10 healthy volunteers (the inclusion criteria of all participants were the same as before) from the Second Hospital Affiliated of Kunming Medical University and anticoagulated them with EDTA to ensure the stability of RNA in the blood samples. Any healthy volunteers had no known disease. Fresh anticoagulated whole blood was added in 3 mL to a 15 mL centrifuge tube containing 3 mL of human peripheral blood lymphocyte isolate (Solarbio, Beijing, China), and the mixed fluid was centrifuged using a high-speed centrifuge (Thermo Fisher Scientific, Waltham, MA, USA) based on 450 g × 30 min at room temperature. After centrifugation, the mixture would be clearly stratified: the top layer is the diluted plasma layer, the middle layer is the clear isolate layer, the white film layer between the plasma and the isolate is the lymphocyte layer, and the bottom of the centrifuge tube is the red blood cells and granulocytes. Subsequently, the cells of the white membrane layer were carefully aspirated into a 15mL clean centrifuge tube and washed with 10mL PBS. The supernatant was discarded and the cells were resuspended in 5 mL of PBS, centrifuged at 250 g for 10 min, and the supernatant was discarded after two repetitions.

RNA Isolation, Reverse Transcription, and Real-Time PCR

Total RNA was isolated from 10 paired samples using TRIzol Reagent (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s instructions. Then the concentration and purity of the RNA solution was quantified using a NanoDrop 2000 nucleic acid protein quantifier (Thermo Fisher Scientific, Waltham, MA, USA). qRT-PCR was performed as described previously.¹⁸ Briefly, the extracted RNA was reverse-transcribed to cDNA using the SureScript-First-strand-cDNA-synthesis-kit (Genecopoeia, Guangzhou, China) prior to qRT-PCR. The qRT-PCR reaction consisted of 4 µL of reverse transcription product, 2 µL of 5X BlazeTaq qPCR Mix (Genecopoeia, Guangzhou, China), 0.5 µL each of forward and reverse primer, and 3-µL nuclease-free water. PCR was performed in a CFX96 Touch instrument (Bio-Rad, Hercules, CA) under the following thermal cycle protocol: (1) 95°C for 30s; (2) 95°C for 10s, 60°C for 20s, 40×; (3) 72° for 30s. The GAPDH gene served as an internal control. RNA levels were calculated for DN samples and normal samples using the 2^−ΔΔCt method. Primers for amplification were purchased from Tsingke (Beijing, China). Primer sequences used for qRT-PCR were shown in Table 1.

Table 1 Primer sequences

Protein Extraction and Western Blot

Western blot analysis was conducted with an SDS-PAGE electrophoresis system (Bio-Rad Laboratories, Hercules, CA, USA) to quantify changes in protein levels. Briefly, total protein was extracted from the PBMC using RIPA buffer containing protease inhibitors and phosphatase inhibitors. Protein samples were separated by SDS-PAGE and transferred onto nitrocellulose membranes, which were subsequently blocked for 12h at 4°C with 5% skimmed milk containing TBST solution. Antibodies against DDX58 (107kDa; ab180675; Abcam), SAMD9L (185kDa; PA5-53994; Thermo), TLR6 (92kDa; ab37072; Abcam), and β-actin (42kDa; ab8226; Abcam) were used as primary antibodies. The samples were incubated with horseradish peroxidase-conjugated secondary antibodies at room temperature (RT) for 1 h. The membrane was imaged by QuickGel 6100 (Monad, Wuhan, China).

Statistical Analysis

The bioinformatic analysis part of this study was performed in the R software. Statistical plots of PCR and Western blot were completed by GraphPad Prism 8.0, and significant P values were calculated by unpaired t-test, where * indicates P < 0.05, ** indicates P < 0.01, and **** indicates P < 0.0001. P < 0.05 was considered statistically significant if not otherwise stated.

Results

Characterization of DE-mRNAs and -lncRNAs Associated with DN

We identified a total of 185 DE-mRNAs between the DN and normal groups, of which 150 were up-regulated genes and the remaining 35 as down-regulated genes (Figure 1A; Supplementary Table 1). Meanwhile, a further 88 DE-lncRNAs were appraised, including 68 up-regulated and 20 down-regulated (Figure 1B; Supplementary Table 2). Additionally, the expression patterns of DE-mRNAs and -lncRNAs in the two groups were graphically visualized with the clustering heatmap (Figures 1C and D). These differentially expressed genes would be available for the subsequent analysis.

Figure 1 Differential expression analysis. (A) volcano plot of 185 DE-mRNAs between the DN and normal groups. (B) volcano plot of 88 DE-lncRNAs between the DN and normal groups. (C and D) clustering heatmap showed the expression pattern of DE-mRNAs (C) and DE-lncRNAs (D) respectively in the two groups (DN and control).

Screening and Functional Enrichment Analysis of Hub Genes

With the miRbase database and miranda software, we constructed a DE-mRNA-miRNA-DE-lncRNA network covering 290 nodes and 496 edges (Figure 2A; Supplementary Table 3). Further, we predicted the protein interactions of 93 DE-mRNAs in the above network utilizing the STRING website and performed visualization in Cytoscape, ultimately, we obtained a PPI network consisting of 42 nodes and 84 edges (Figure 2B; Supplementary Table 4). Subsequently, 10 genes defined as hub genes were filtered from the above PPI network according to the CytoHubba plug-in (Supplementary Table 5) and displayed in Figure 2C, which had a total of 37 edges.

Figure 2 Establishment of ceRNA network and screening of the hub genes. (A) the DE-mRNA-miRNA-DE-lncRNA network (ceRNA) covers 290 nodes and 496 edges. (B) the PPI network consists of 42 nodes and 84 edges, utilizing the STRING website to predict the protein interactions of 93 DE-mRNAs and performed visualization in Cytoscape. (C) 10 hub genes were filtered from the PPI network. (D and E) GO (D) and KEGG (E) analysis of the 10 hub genes.

To reveal the potential functions of hub genes in ND, GO and KEGG analysis were imperative (Figures 2D and E). We found that hub genes appeared significantly enriched in terms associated with the response to viral infection. For example, the ‘defense response to virus’, ‘response to virus’, “cytoplasmic pattern recognition receptor signaling pathway in response to virus”, “cellular response to virus”, “negative regulation of viral process”, “regulation of defense response to virus by host”, “negative regulation of viral genome replication”, “viral genome replication”, etc. This seemed to coincide with previous reports that viral infections (eg, EBV, HCV, HBV, etc.) induce the development of DN.^19–21 Similarly, KEGG also suggested that the hub genes were tightly coupled with “Epstein-Barr virus infection”, “Hepatitis C”, and “Hepatitis B”. Intendedly, hub genes were implicated in a variety of processes related to kidney development, including “metanephric mesenchyme development”, “nephron tubule epithelial cell differentiation”, “negative regulation of kidney development”, “renal vesicle morphogenesis”, and “metanephric nephron tubule development”. Furthermore, the hub gene may also influence the progression of DN by regulating the immune response (“regulation of innate immune response”, “regulation of immune effector process”, “activation of innate immune response”, “positive regulation of innate immune response”, “positive regulation of immune effector process”, etc.) and the status of immune cells (“regulation of dendritic cell cytokine production”, “regulation of granulocyte macrophage colony-stimulating factor production”, “positive regulation of macrophage activation”, etc.). Otherwise, “regulation of inflammatory response”, “NLRP3 inflammasome complex assembly”, “regulation of NLRP3 inflammasome complex assembly”, and “positive regulation of inflammatory response” being also remarkably enriched. The exhaustive results of the GO and KEGG analyses were available in Supplementary Tables 6 and 7, respectively. These results suggested that the hub genes may play a role in the pathogenesis of DN by regulating processes related to viral infection, kidney development, and inflammatory and immune responses.

Diagnostic Biomarkers Filtering for DN from Hub Genes

The minimum lambda value was visible in 3 following the cvfit plot, which indicated that the 3-hub gene-based diagnostic model was optimized in terms of both accuracy and simplicity (Figure 3A). Therefore, only DDX58, SAMD9L, and TLR6 were included (Figure 3B). A ROC curve analysis was conducted to evaluate the sensitivity and specificity of the 3-hub gene-based diagnostic model for differentiating DN from the normal group (Figure 3C). By calculating the AUC, it was revealed that the 3-hub gene signature could accurately classify DN patients from normal individuals (AUC = 1). Also, individual ROC curve analysis of the three hub genes mentioned above indicated that these genes could also be independently distinguished between DN patients and normal subjects (AUC = 1 for all; Figure 3D). Hence, these genes were considered as diagnostic biomarkers for DN and served for subsequent analysis.

Figure 3 LASSO regression analysis of the 10 hub genes. (A and B) LASSO regression analysis screened the 3-hub gene-based diagnostic model. (C) ROC curve analysis showed the sensitivity and specificity of the 3-hub gene-based diagnostic model for differentiating DN from the normal group. (D) individual ROC curve analysis of the 3 hub genes.

Functional Enrichment Analysis of Diagnostic Biomarkers

To further identify the potential pathways involved in diagnostic biomarkers, we performed single-gene GSEA-KEGG on them and the results were displayed in Figures 4A–C. Collectively, all three diagnostic biomarkers were associated with viral infection (“Herpes simplex virus 1 infection”, “Hepatitis C”, “Epstein-Barr virus infection”, etc.). Surprisingly, we found that the “NOD-like receptor signaling pathway”, which is involved in the immune regulation of DN,^22,23 was significantly enriched in the above results. Notably, SAMD9L was involved in “Diabetic cardiomyopathy”. Similarly, DDX58 had been linked to cardiac disease pathways (“Dilated cardiomyopathy” and “Arrhythmogenic right ventricular cardiomyopathy”). Detailed enrichment results for each diagnostic biomarker were reported in Supplementary Tables 8–10. In addition, we constructed a ceRNA network encompassing 16 KEGG pathways based on three diagnostic biomarkers, which included three diagnostic biomarkers, three miRNAs, and seven DE-lncRNAs (Figure 4D).

Figure 4 Functional enrichment analysis of the 3 hub genes. (A–C) single-gene GSEA-KEGG was performed on the 3 hub genes, DDX58 (A), SAMD9L (B), and TLR6 (C and D) a network containing four parts, 3 hub genes and 16 KEGG pathways based on those genes, 3 miRNAs, and 7 DE-lncRNAs.

Expression of Diagnostic Markers Was Frequently Elevated in DN Samples

To verify the expression of three diagnostic biomarkers (DDX58, SAMD9L, and TLR6) in DN, qRT-PCR was used to detect the mRNA levels of DDX58, SAMD9L, and TLR6 in PBMC obtained from 10 cases of DN and healthy subjects. The results showed that the mRNA levels of DDX58, SAMD9L, and TLR6 were significantly higher in the DN group compared with healthy subjects (P < 0.05, Figure 5A–C). Moreover, Western blot showed that protein levels of DDX58, SAMD9L, and TLR6 were significantly increased in DN samples compared to healthy subjects (Figure 5D). In summary, our results showed that DDX58, SAMD9L, and TLR6 were upregulated in DN patients at both transcriptome and protein levels compared to healthy controls.

Figure 5 Experimental validations of the 3 hub genes. (A–C) qRT-PCR showed the mRNA levels of DDX58 (A), SAMD9L (B), and TLR6 (C and D) Western blot showed the protein levels of the 3 genes. N represented the control group and D represented the DN group. * P<0.05 versus control, **P<0.01 versus control, ****P<0.001 versus control.

Discussion

As we all know, DN is not only the major cause of ESKD but also a serious complication of diabetes. Early DN may not have any clinical symptoms, but only glomerular hyperfiltration.²⁴ Although renal biopsy is the “gold standard” for the diagnosis of renal disease, glomerular changes in the early stage of diabetes are non-specific and renal biopsy as an invasive operation does have lots of risks, thus it may not be widely accepted in this period.¹⁴ In addition, there is no other test to replace renal biopsy at present. However, as for DN, diagnosis and cure earlier are of great significance. Drug intervention in the early stage could delay the progression of renal failure²⁵ and prevent the occurrence of other adverse events such as cardiovascular disease²⁶ (CVD), thereby improving patients’ quality of life. Therefore, it is urgent to find a new biomarker of early DN.

In this study, 185 DE-mRNAs and 88 DE-lncRNAs were obtained by sequencing the samples of DN and healthy volunteers. Because of the large number of differential expression genes, it is very hard to study these genes one by one. We all know that lncRNA is a kind of non-coding RNA.²⁷ Therefore, the study of lncRNA function is usually carried out in lncRNA-RNA, lncRNA-DNA, lncRNA-protein.²⁸ In recent years, many studies have suggested that miRNA could participate in the regulation of many important cellular activities, such as cell proliferation, differentiation, epithelial-mesenchymal transformation and so on,^29,30 making miRNA gradually has been a hot area of research. And then, it was found that lncRNA could serve as a miRNA sponge, reducing the degradation or transcriptional inhibition of mRNA by binding to miRNA, which is often referred to as ceRNA.³¹ Therefore, we hope to further screen the genes which are more relevant to disease among those differential expression genes we got by constructing such a sort of ceRNA network (lncRNA-miRNA-mRNA). Fortunately, we linked the ceRNA network with the PPI network, and 10 hub genes were screened, then the functional enrichment analysis was performed. The diagnostic model constructed by LASSO regression analysis was used to gain the three characteristic genes (DDX58, SAMD9L, and TLR6), which are most related to DN. Finally, the single-gene GSEA analysis of the three genes was carried out, and the gene-miRNA-lncRNA-KEGG pathway was constructed.

We showed great interest in the three genes we finally got, the reason was that these genes showed great sensitivity and specificity to distinguish between DN and non-DN. DDX58, called DExD/H-box helicase 58, is deemed to take part in viral ds RNA recognition and the regulation of the antiviral innate immune response. SAMD9L, called sterile alpha motif domain containing 9 like, is regarded as playing a key role in cell proliferation and the innate immune response to viral infection. TLR6, called toll-like receptor 6, is considered to play a fundamental role in pathogen recognition and activation of innate immunity. In consequence, we speculated that these three genes maybe new biomarkers of DN.

As mentioned before, by observing and comparing the single-gene GSEA results, we found that, without exception, all the genes were associated with “viral infection” and coincide with a “NOD-like receptor signaling pathway”.

It’s well known that viral infection is one of the causes of many kidney diseases. From the latest reports, viruses related to kidney diseases include hepatitis B virus,³² hepatitis C virus,³³ human immunodeficiency virus,³⁴ cytomegalovirus, Epstein-Barr virus, parvovirus³⁵ and so on. Covid-19, which is now prevalent all over the world, has also been found to cause kidney damage.³⁶ The pathogenesis of virus infection-associated glomerular disease is still unclear. Some studies suggested that it may be due to virus infection resulting in renal tissue damage, or through viral antigens and corresponding antibodies to form the circulating immune complexes (CIC) deposition, and then activated the specific immunity mechanism.³⁷ Consequently, we will discuss the aspect of immune defense principally.

Some studies have shown that virus infection was the main cause of autoimmunity.³⁸ The innate immune system’s recognition of microbial pathogens was mediated by pattern-recognition receptors (PRR), while PRR detected pathogen-associated molecular patterns (PAMPs) or danger-associated molecular patterns (DAMPs), such as viral nucleic acids.³⁹ At present, the widely accepted PRR mainly includes Toll-like receptors (TLRs), C-type lectin receptors (CLRs), RIG-I-like receptors (RLRs), Nod-like receptors (NLRs), and AIM-2 like receptors.⁴⁰ It’s also suggested that DDX58, SAMD9L, and TLR6 could indeed participate in innate immune responses through the NOD-like receptor signaling pathway.

NLRs are a cluster of multidomain proteins containing an alterable amino-terminal (N-terminal) domain, a central NOD module, and carboxy-terminal leucine-rich repeats (LRRs).⁴¹ There are two prototypic NLRs, NOD1 and NOD2. On the one hand, they could sense the cytosolic presence of the bacterial peptidoglycan fragments which escaped from endosomal compartments, thus driving the activation of nuclear factor kappa B (NF-κB) and interferon regulatory factor (IRF), cytokine production, and apoptosis. On the other hand, they could also induce caspase-1 activation, to regulate the maturation of IL-1B, and IL-18 and drive pyroptosis. Some studies have indicated that the activation of NLRs will lead to the formation of inflammasome.⁴² However, there’re a large number of studies published regarding the effect and mechanism of NLRP inflammasomes in the pathogenesis and progression of DN,^43,44 so it needs no further elaboration here.

The following limitations exist in this study. First, a single source of our research data and the limited sample size could cause some deviations. Secondly, the study needs to be further supplemented to reveal the mechanism of the key genes. Lastly, the diagnostic model remains unproven in more clinical samples.

In general, we obtained 10 hub genes by connecting the ceRNA network with the PPI network. Then the diagnostic model was constructed by LASSO regression analysis to screen three characteristic genes related to the disease. But this was only our guess, so we further verified the genes by regaining other clinical samples of DN, and confirmed that these three genes were indeed highly expressed in the DN, while it was opposite in the control. This obviously differential expression phenomenon of the disease maybe indicate that these three genes are closely related to the occurrence of DN. Through the enrichment analysis of these genes, the key signaling pathway was obtained. We guessed that, on the one hand, these three genes may be able to identify renal damage caused by a viral infection, and the other side may be able to participate in the progression of DN through NOD-like receptor signaling pathway-mediated innate immune responses and inflammatory. However, the more specific mechanism remains to be further studied.