Peripheral Blood circRNA Microarray Profiling Identities hsa_circ_0001831 and hsa_circ_0000867 as Two Novel circRNA Biomarkers for Early Type 2 Diabetic Nephropathy

Introduction

Diabetic nephropathy (DN) occurs in approximately 40% of Type 2 diabetes mellitus (T2DM),¹ with the progress of T2DM; early type 2 diabetic nephropathy (ET2DN) eventually develops into end-stage renal disease.^2,3 Thus, it is particularly crucial to enhance the prevention and treatment of ET2DN and discover valid biological diagnostic biomarkers.

Circular RNAs (circRNAs) are a group of newly discovered non-coding RNAs.⁴ In recent years, there has been accumulating evidence that circRNAs can regulate host gene’s splicing and transcription and act as microRNAs and protein sponges, as well as moderate protein translation and protein functions.⁵ CircRNAs have been identified to play an active role in treating a broad range of diseases.^6–8

The pathological manifestations associated with DN include glomerular thylakoid cells hypertrophy and hyperplasia, extracellular matrix (ECM) accumulation, podocytes function impairment and renal fibrosis.^9,10 Some studies have been partially demonstrated that circRNAs can act as microRNA sponges to influence the pathogenesis of DN, for instance, An et al have shown that hsa_circ_0003928 can bind to miR-151-3p and moderate Annexin A2 (ANXA2) expression via the microRNAs sponge, thereby attenuating high glucose-induced apoptosis in podocytes and glomerular thylakoid cells.¹¹ circRNA_0080425 is able to target miR-24-3p to regulate fibroblast growth factor 11 (FGF11) expression levels, which led to fibrosis and promoted DN progression.¹² miR-24-3p play an important role in regulating cell proliferation, apoptosis and fibrosis.^13–15 A study in mice with DN showed that circRNA_15698, which is highly expressed in DN mice, functions as a sponge for miR-185 and upregulates the expression of TGFB1, thereby promoting the production of ECM proteins.⁹

In addition, some circRNAs have been reported to be particular in DN via inflammatory pathways, for example, circRNA_ITCH improved the renal inflammation in streptozotocin-induced diabetic mice,¹⁶ knockdown of circLRP6 decreased high glucose-induced NF-κB signaling pathway in mesangial cells,¹⁷ circRNA_WBSCR17 accelerated inflammatory responses in high glucose-induced kidney tubular cells.¹⁸

There has been a study on the expression profile of circRNAs using kidney tissue from DN mice, and differential expression circRNAs (DECs) may be candidated as some biomarkers for DN.¹⁹ However, there has been no study on circRNAs as diagnostic biomarkers for ET2DN in humans yet. The aim of this study was to find new diagnostic biomarkers for ET2DN by profiling the DECs in peripheral blood from patients with ET2DN.

In this study, we performed circRNA microarray assay on peripheral blood samples from patients with ET2DN and T2DM without DN and age-sex matched normal controls to screen for DECs. Partial circRNAs were validated by quantitative real‐time reverse transcription‐PCR (qRT‐PCR), and further bioinformatic analysis was performed to explore their biological functions. ROC curves were made to evaluate clinical diagnostic value of DECs.

Subjects and Methods

Study Subjects

A total of 60 subjects, including 20 ET2DN patients, 20 age- and sex-matched T2DM without DN patients and 20 age- and sex-matched normal control subjects, were enrolled in this study. All ET2DN and T2DM groups were presented in the Department of Endocrinology of the Third Xiangya Hospital of Central South University from November 2021 to June 2022, and all control subjects volunteered to participate in this study.

The ET2DN inclusion criteria included patients who were firstly tested by urinary albumin excretion rate (UAER) in 20–200 μg/min and were further verified by UAER in 30–300 mg/24h for twice in different days. Patients in ET2DN and T2DM group were confirmed according to World Health Organization criteria.²⁰ All subjects with serious infectious diseases, autoimmune diseases, hematological diseases, malignant diseases, fever were excluded, and control subjects had no history of diabetes mellitus or renal disease. The clinical characteristics of the participants included age, gender, duration of diabetes, systolic blood pressure (SBP), diastolic blood pressure (DBP), height, weight, BMI, waist circumference (WC), hip circumference (HC). In addition, biochemical parameters including fasting blood glucose (FBG), HbA1c, creatinine (Cr), UAER, total cholesterol (TC), low-density lipoprotein cholesterol (LDL-C), high-density lipoprotein cholesterol (HDL-C), triglycerides (TG) and haemoglobin (Hb) levels were measured.

This study was approved by the Ethics Committee of The Third Xiangya Hospital of Central South University, Changsha, China, and was performed according to the guidelines outlined in the Declaration of Helsinki. All participants provided written informed consent form.

Peripheral Blood Samples Collection

Two milliliters of peripheral venous blood samples from different groups’ patients were collected. The samples were then randomly selected from 3 ET2DN, 3 age- and sex-matched T2DM and 3 age- and sex-matched control subjects for RNA extraction and circRNA microarray analysis.

CircRNA Microarray Analysis

Total RNA was isolated from peripheral blood samples with the TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The NanoDrop ND-1000 (Nanodrop Technologies, Wilmington, DE) was used to quantify the purity and concentration of total RNA samples. Preparation of samples and microarray hybridization was conducted following Arraystar’s standard protocol. Briefly, digestion of total RNA with RNase R (Epicentre, Madison, WI, USA) was performed to remove linear RNA and enriched circRNAs. The abundant circRNA is in turn amplified and transcribed into fluorescent cRNA by the random initiation method (Arraystar Super RNA Label Kit; Arraystar), and tagged cRNAs were hybridized to Arraystar Human CircRNA arrays (8x15K, Arraystar). Once the slides were cleaned, the arrays were scanned using an Agilent G2505C scanner (Agilent Technology, Santa Clara, CA, USA).

Then, we featured extraction software from Agilent (version 11.0.1.1) for analysis of the acquired array images. The R software limma-package was used for the quantitative normalization and subsequent data processing. CircRNAs with significant differential expression between the three groups were shown by Volcano Plot. The DECs between two samples were identified by fold-change filtering. Hierarchical clustering showed distinguishable expression patterns of circRNAs between samples. The DECs between two groups were defined as | log2 fold change| ≥ 1 and p-value < 0.05.

Gene ontology (GO) analysis of DECs according to their host gene was carried out by GO seq R package, where the ontology covered three areas: biological processes (BP), cellular components (CC) and molecular functions (MF). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of relevant biological pathways was performed by KEGG database (http://www.geneontology.org). The p-value produced by top GO denotes the significance of GO term enrichment in differential circRNAs.

qRT-PCR Validation of Candidate circRNAs

The DECs between the each groups were validated by qRT‐PCR in accordance with the raw data analysis. Total RNA was reversed transcribed into cDNA using Revert Aid TMM-MuLV reverse transcriptase (Invitrogen, Carlsbad, CA, USA) with random primers according to the manufacturer’s instructions. RT-PCR was performed at real time on a VIIa 7 real-time polymerase chain reaction system (Applied Biosystems, New York, NY, USA). The thermal cycling conditions for the PCR reaction were 95 °C for 10 min, followed by 45 cycles of 95°C for 10s, 60°C for 60s and 95 °C for 15s. The GAPDH was used as the internal control. The primer list is provided in Supplementary Table 1. The relative expression level was calculated using the 2^−∆∆ct method.

The Prediction of circRNA–miRNA–mRNA Networks

To investigate the biological potential function of the target circRNAs, the bioinformation method was used for the prediction of circRNA–microRNA–mRNA interaction network. Firstly, the circRNA–microRNA interaction was predicted with Arraystar’s home-made miRNA target prediction software based on TargetScan (https://www.targetscan.org/) and miRanda (http://www.microrna.org/microrna/home.do). DIANA TOOLS (http://diana.imis.athena-innovation.gr/DianaTools/index.php) software was used to predict the microRNA and mRNA interaction; the circRNA–microRNA–mRNA network was displayed using Cytoscape 2.8.2.²¹

Receiver Operating Characteristic (ROC) Curves

ROC curves were plotted using the ET2DN group as the state variable and the relative expression level of hsa_circ_0001831 and hsa_circ_0000867 as the test variable to assess the diagnostic performance of DECs against ET2DN. The ROC curves and the area under the curve (AUC) were analyzed, and sensitivity as well as specificity was calculated using 95% confidence intervals (Cis), with AUC > 0.5 representing diagnostic value.

Statistics Analysis

Statistical analyses were performed using SPSS 26.0 (SPSS, Inc., IL, USA) and GraphPad Prism 8.0 (GraphPad, Inc., CA, USA). Clinical data and results from experiments were presented as means ± standard deviation or numbers with percentages. Shapiro Wilk test was used to verify whether random sample data conformed to normal distribution. Normally distributed data between groups were calculated by independent sample t-test and one-way analysis of variance (ANOVA), non-normally distributed data were tested using Mann–Whitney U-test. P ≤ 0.05 was considered as statistically significant difference. qRT-PCR results were compared using the unpaired t‐test. Multiple logistic regression was used to analyze the correlation between circRNAs and clinical indicators.

Results

Characteristics of the Study Population

The basic clinical characteristics of the subjects are revealed in Table 1. As shown, the duration of diabetes, SBP, DBP, FBG, HbA1c and HDL-C levels had significant difference in ET2DN vs control group; the ET2DN group had substantially higher UAER than the other two groups (P = 0.014). However, there were no statistical differences in age, height, weight, WC, HC, BMI, Cr, TC, TG, LDL-C and Hb levels across the three groups.

Table 1 Characteristics and Clinical Features of Study Population

DECs Between Groups

High-throughput human circRNA microarray revealed the circRNA expression profiles in the peripheral blood of three groups’ subjects. By quantile normalization of raw data for further analyses, the Box Plot showed that the distribution of log 2 ratios was similar in all the samples (Figure 1A). Hierarchical clustering analysis showed that the circRNA expression pattern was significantly different between the three groups (Figure 1B). Fold-change filtering of DECs was showed in scatter plot (Figure 1C-E). In addition, Volcano Plots was used to visualize the variation in circRNA expression with statistical significance between the three groups (fold change >1.0 and P < 0.05) (Figure 2A-C).

Figure 1 DECs in ET2DN, T2DM and control groups. (A) The Box Plot displayed the normalized intensity distribution of all data. (B) Hierarchical clustering for DECs between ET2DN, T2DM, and control. High relative expression in red, low relative expression in green. Scatter plots were performed to analyze DECs in ET2DN vs T2DM group (C), T2DM vs control group (D) and ET2DN vs control group (E). CircRNAs above the top green line and below the bottom green line indicated more than two-fold change in circRNAs between the two groups.

Figure 2 Volcano Plots in ET2DN, T2DM and control groups indicated the DECs between ET2DN vs T2DM group (A), T2DM vs control group (B) and ET2DN vs control group(C). Fold change values are marked by the vertical lines. A green vertical line corresponds to a 1-fold change, while a horizontal green line marks a p-value <0.05, red point represents DECs.

According to the circRNA microarray results, a total of 586 circRNAs were differentially expressed between the ET2DN and T2DM group, of which 249 were upregulated and 337 were downregulated. Top 10 upregulated and downregulated circRNAs were shown in Tables 2 and 3, among them, hsa_circ_0001831 and hsa_circ_0044235 were the most upregulated and downregulated circRNAs between the ET2DN and T2DM group, respectively. In addition, 176 circRNAs were upregulated and 101 circRNAs were downregulated in T2DM vs control group; top 10 upregulated and downregulated circRNAs are shown in Supplementary Table 2 and Supplementary Table 3. Fifty-seven circRNAs were upregulated and 5 circRNAs were downregulated in ET2DN vs control group; top 10 upregulated and 5 downregulated circRNAs are shown in Supplementary Table 4 and Supplementary Table 5, respectively.

Table 2 Upregulated of circRNAs in ET2DN Vs T2DM Group

Table 3 Downregulated of circRNAs in ET2DN Vs T2DM Group

We found that majority of the DECs of three groups on human chromosome distribution were mainly located on chromosome 1 (Supplementary Figure 1A, B, C).

GO and KEGG Pathway Analysis for Host Genes of circRNAs

To predict the function enrichment of DECs, GO and KEGG analyses were performed. GO analysis of upregulated circRNAs in ET2DN vs T2DM group showed that the DECs were enriched in small GTPase mediated signal transduction and positive regulation of GTPase activity and so on (Figure 3A). The analysis of downregulated circRNAs in ET2DN vs T2DM group suggested that DECs were more relevant to biological processes such as symbiont process and viral process (Figure 3B). In addition, GO analysis of T2DM vs control group was performed (Supplementary Figure 2A, B); however, GO analysis could not be enriched because too few DECs in ET2DN vs control group. KEGG analysis indicated that the upregulation of circRNAs remarkably enriched pathways including morphine addiction, human cytomegalovirus infection, and downregulation mainly focused on RNA transport, mTOR signaling pathway (Figure 3C and D). Pathway analysis of T2DM vs control group was analyzed (Supplementary Figure 3A, B). As mentioned above, pathway analysis cannot be enriched because there are too few DECs in ET2DN vs control group.

Figure 3 GO and pathway analysis of DECs for their host genes in ET2DN vs T2DM group. GO analysis of upregulated (A) and downregulated (B) DECs in ET2DN vs T2DM group. Pathway analysis of upregulated (C) and downregulated (D) DECs in ET2DN vs T2DM group. (E) The common circRNAs were analyzed in ET2DN vs T2DM group and ET2DN vs control group.

To explore whether special circRNAs were upregulated or downregulated in ET2DN compared to T2DM and control group, the common circRNAs from ET2DN vs T2DM group and ET2DN vs control group were further analyzed. Total 3 circRNAs, named hsa_circ_0000867, hsa_circ_0022430 and hsa_circ_0063812, were upregulated in ET2DN compare to T2DM and control group; however, hsa_circ_0009035 was downregulated (Figure 3E).

Validation of DECs by qRT-PCR

Based on the results of circRNA microarray, six candidate circRNAs were selected for validation by qRT-PCR, including the highest expression (hsa_circ_0001831) and the lowest expression (hsa_circ_0044235 and hsa_circ_0061260) in ET2DN vs T2DM group; the highest expression (hsa_circ_0044235) and the lowest expression (hsa_circ_0001831) in T2DM vs control group; the highest expression (hsa_circ_0000658) and the lowest expression (hsa_circ_0063503) in ET2DN vs control group; the highest expression of circRNAs in ET2DN vs T2DM and control group (hsa_circ_0000867). The Results of RT‑qPCR suggested that hsa_circ_0001831 was significantly upregulated in ET2DN vs T2DM group. In addition, although hsa_circ_0001831 did not be shown in ET2DN vs control group in our circRNA microarray data, our RT-PCR results showed that hsa_circ_0001831 was also upregulated in ET2DN vs control group (Figure 4A). In common expression data, hsa_circ_0000867 was verified as higher expression in ET2DN vs T2DM and control group (Figure 4B). However, hsa_circ_0061260 and hsa_circ_0044235 did not have difference in ET2DN vs T2DM group (Figure 4C and D); hsa_circ_0063503 and hsa_circ_0000658 also had no statistically significant difference between the ET2DN vs control group (Figure 4E and F); hsa_circ_0044235 had no statistically significant difference in T2DM vs control group (Figure 4G). Altogether, consistent with the microarray results, the expression level of hsa_circ_0001831 and hsa_circ_0000867 may be two diagnosis biomarkers in ET2DN in peripheral blood. Therefore, hsa_circ_0001831 and hsa_circ_0000867 were selected for further study.

Figure 4 Verification of DECs by qRT-PCR. The expression of hsa_circ_0001831(A) and hsa_circ_0000867(B) in three groups; The expression of hsa_circ_0061260(C) and hsa_circ_0044235(D) in ET2DN vs T2DM group; The expression of hsa_circ_0063503(E) and hsa_circ_0000658 (F) in ET2DN vs control group; The expression of hsa_circ_0044235 (G) in T2DM vs control group. n=20, *p<0.05, **p<0.01, ***p<0.001.

Prediction of the circRNA–microRNA–mRNA Interaction Network

To further clarify the function of hsa_circ_0001831 and hsa_circ_0000867, TargetScan and miRanda were used to predict circRNA–microRNA interaction network, and DIANA TOOLS software was used to predict the microRNA–mRNA interaction. This indicated that hsa-miR-3151-5p, hsa-miR-6791-5p, hsa-miR-939-5p, hsa-miR-637 and hsa-miR-7974 were predicted as target microRNAs of hsa_circ_0001831, and has-miR-146a-3p and has-miR-877-3p were predicted as target microRNAs of hsa_circ_0000867. The interaction network of hsa_circ_0001831 and hsa_circ_0000867 with microRNA–mRNA is shown in Figure 5A and B.

Figure 5 Prediction of hsa_circ_0001831 and hsa_circ_0000867–microRNA–mRNA interaction network. The possible microRNAs of hsa_circ_0001831 (A) and hsa_circ_0000867 (B) were predicted by Arraystar’s home-made software based on TargetScan and miRanda, next, DIANA TOOLS software was used to predict target mRNAs of microRNA, Cytoscape 2.8.2 was used to display circRNA–microRNA–mRNA interaction network.

The Value of hsa_circ_0001831 and hsa_circ_0000867 in the Diagnosis of ET2DN

To determine the diagnostic value of hsa_circ_0001831 and hsa_circ_0000867 in ET2DN, ROC curves were drawn and analyzed. The results showed that hsa_circ_0001831 had a sensitivity of 0.85, a specificity of 0.85, an AUC value of 0.95 and a maximum Jordan index of 0.75 for the diagnosis of ET2DN (Figure 6A). In addition, hsa_circ_0000867 had a sensitivity of 0.6, a specificity of 0.95, an AUC value of 0.95 and a maximum Jordan index of 0.55 for the diagnosis of ET2DN (Figure 6B). Based on the Jordan index, the optimal exploration cut-off point for hsa_circ_0001831 was 1.55 and hsa_circ_0000867 was 1.14, respectively. Together, hsa_circ_0001831 and hsa_circ_0000867 showed high diagnostic value for ET2DN.

Figure 6 ROC curve analysis of diagnostic value in hsa_circ_0001831 (A) and hsa_circ_0000867 (B) in ET2DN patients.

Discussion

DN is one of the chronic complications of diabetes, which causes end-stage renal disease.^22,23 Some studies have revealed that circRNAs can be used for diagnostic biomarkers for some diseases.^24–26 CircRNAs may be involved in the development of DN and could be a new therapeutic target for DN, for instance, circ_0000285 caused podocyte injury by sponging miR-654-3p and activating MAPK6 in DN.²⁷ However, there was no study on the expression profile of circRNAs in ET2DN. Therefore, our aim was to explore the expression profile of circRNAs in patients’ peripheral blood with ET2DN, and we tried to find out some new diagnostic biomarkers for ET2DN.

In the present study, circRNA microarray was used for detecting the DECs in peripheral blood of patients with ET2DN. Our results indicate that there were 586 DECs in ET2DN vs T2DM group and 62 DECs in ET2DN vs control group. These results implied that the DECs in ET2DN patients were significantly different from T2DM and control group, and these circRNAs might be strongly associated with ET2DN.

Next, we performed GO and KEGG analysis on the DECs to understand the possible functions of circRNAs in ET2DN. The results suggested that the DECs in ET2DN were closely associated with multiple biological processes such as GTPases and protein modifications, and the pathways were enriched in some biological pathways like RNA transport, mTOR signaling pathway. Several pathways have been reported to be related to DN, including MAPK and mTOR signaling pathway.²⁸ We speculated that some circRNAs might be involved in the modification of ET2DN through these pathways.

To further validate the DECs, six candidate DECs were detected by qRT-PCR. The expression of hsa_circ_0001831 and hsa_circ_0000867 were significantly upregulated in ET2DN compared to both of T2DM and control group. Therefore, hsa_circ_0001831 and hsa_circ_0000867 may be two diagnostic biomarkers in ET2DN. hsa_circ_0001831 was produced by its host gene, scribble planar cell polarity protein (SCRIB), which is a tumor suppressor, hsa_circ_0001831 could decrease SCRIB mRNA expression but did not affect mRNA splicing, in addition, hsa_circ_0001831 was upregulated in breast cancer and promoted breast cancer cell proliferation, migration and invasion.²⁹ Moreover, SCRIB, as an essential gene, was upregulated in podocytes according to the Human Protein Atlas (HPA) database (https://www.proteinatlas.org).³⁰

However, there were no study about hsa_circ_0000867, but its host gene, ornithine decarboxylase antizyme 1 (OAZ1) was significantly associated with T2DM duration and correlated with HbA1c.^31,32

CircRNAs indirectly promoted target mRNA expression through binding to microRNAs and repress microRNAs expression.³³ Therefore, we predicted the hsa_circ_0001831 and hsa_circ_0000867-microRNA–mRNA interaction network. In predicted target microRNAs of hsa_circ_0000867, hsa-miR-877-3p were upregulated in urinary exosomal of T2DN.³⁴ In addition, hsa-miR-146a-3p was strongly upregulated in DN patients compared to control group;³⁵ however, another study showed that inhibition of hsa-miR-146a increases podocyte injury in diabetic glomerulopathy³⁶ and inhibited inflammatory and reactive oxygen species in mice DN.^37,38 Moreover, inflammatory pathways play a critical role in T2DN, except for hsa-miR-146a, some predicted microRNAs were related to inflammatory process in other diseases, for instance, hsa-miR-939 negatively regulated pro-inflammatory genes in human monocytic leukaemia THP-1 cells.³⁹ Furthermore, miR-939-5p was able to attenuate PM2.5-induced endothelial injury via NF-κB p65 signaling pathway in human aortic endothelial cells,⁴⁰ and hsa-miR-637 suppressed vascular inflammation in hypertension through inflammatory mediators.⁴¹ Besides, miR-877-3p negatively regulated inflammatory factor including IL-1β in mesangial cells.⁴² Together, hsa_circ_0001831 and hsa_circ_0000867 might induce ET2DN through these predicted microRNAs. However, more research is needed on the mechanisms by which these DECs affect T2DN.

An assessment of the potential diagnostic value of hsa_circ_0001831 and hsa_circ_0000867 for ET2DN was made. Analysis using ROC curves based on disease diagnosis showed that hsa_circ_0001831 and hsa_had a larger AUC value, better sensitivity and strength, as well as higher predictive accuracy. Therefore, hsa_circ_0001831 and hsa_circ_0000867 could be the two potential new biomarkers for ET2DN screening and diagnosis.

There are some limitations in our study. Firstly, the sample size of this study was relatively small. Secondly, only a few circRNAs from the DECs were validated. Thirdly, the specific mechanism of these DECs in ET2DN remains to be further investigated.

Conclusion

In summary, our study revealed the circRNA expression profiles of ET2DN patients utilizing circRNA microarrays indicating that hsa_circ_0001831 and hsa_circ_0000867 may be used as two biomarkers of ET2DN, which provides a new idea for the diagnosis of ET2DN.