{'fraction_non_alphanumeric': 0.03245829050367676, 'fraction_numerical': 0.02383964576579426, 'mean_word_length': 5.115812379110252, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 4, 'https://': 2, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Urinary sediment CCL5 messenger RNA as a potential prognostic biomarker of diabetic nephropathy

ABSTRACT Background Urinary sediment messenger RNAs (mRNAs) have been shown as novel biomarkers of kidney disease. We aimed to identify targeted urinary mRNAs in diabetic nephropathy (DN) based on bioinformatics analysis and clinical validation. Methods Microarray studies of DN were searched in the GEO database and Nephroseq platform. Gene modules negatively correlated with estimated glomerular filtration rate (eGFR) were identified by informatics methods. Hub genes were screened within the selected modules. In validation cohorts, a quantitative polymerase chain reaction assay was used to compare the expression levels of candidate mRNAs. Patients with renal biopsy–confirmed DN were then followed up for a median time of 21 months. End-stage renal disease (ESRD) was defined as the primary endpoint. Multivariate Cox proportional hazards regression was developed to evaluate the prognostic values of candidate mRNAs. Results Bioinformatics analysis revealed four chemokines (CCL5, CXCL1, CXLC6 and CXCL12) as candidate mRNAs negatively correlated with eGFR, of which CCL5 and CXCL1 mRNA levels were upregulated in the urinary sediment of patients with DN. In addition, urinary sediment mRNA of CXCL1 was negatively correlated with eGFR (r = −0.2275, P = 0.0301) and CCL5 level was negatively correlated with eGFR (r = −0.4388, P < 0.0001) and positively correlated with urinary albumin:creatinine ratio (r = 0.2693, P = 0.0098); also, CCL5 and CXCL1 were upregulated in patients with severe renal interstitial fibrosis. Urinary sediment CCL5 mRNA was an independent predictor of ESRD [hazard ratio 1.350 (95% confidence interval 1.045–1.745)]. Conclusions Urinary sediment CCL5 and CXCL1 mRNAs were upregulated in DN patients and associated with a decline in renal function and degree of renal interstitial fibrosis. Urinary sediment CCL5 mRNA could be used as a potential prognostic biomarker of DN.


INTRODUCTION
The worldwide prevalence of diabetes mellitus (DM) has increased dramatically in the past few decades. Diabetic nephropathy (DN) is an established chronic complication of DM, which occurs in ∼30% of patients with type 1 DM (T1DM) and 40% of type 2 DM (T2DM) patients [1]. The overall prognosis of DN is relatively poor. About 20-30% of individuals with DN will progress to end-stage renal disease (ESRD) even when on optimal glycemic control and therapy with a renin-angiotensinaldosterone system (RAAS) blocker [2]. Estimated glomerular filtration rate (eGFR), an indicator of renal function, has been regarded as the preferred diagnostic and prognostic biomarker of DN in clinical practice or trials [3]. However, the abnormality of eGFR only reflects an advanced stage rather than the early functional changes. Albuminuria also predicts the progression of DN, but it lacks specificity and sensitivity for progressive decline in eGFR and occurrence of ESRD. A proportion of renal impairment happened even before the appearance of albuminuria [4]. The biomarkers currently in use are not accurate enough to predict the risk for the above-mentioned limitations. Thus, it is still a challenge to identify DN patients with a high risk of progression.
Previously, various types of cells have been described in the urinary sediment of patients with renal disease. Quantification of targeted messenger RNAs (mRNAs) in urinary sediments represents a potential non-invasive diagnostic method for multiple types of kidney injuries. The detection of podocyte markers [5], molecules associated with epithelial-mesenchymal transition [6] and inflammatory cytokines [7,8] has emerged as a potential approach to identify novel biomarkers for diagnosis of DN. In recent years, tools of bioinformatics have also been applied to screen urinary mRNAs in DN, providing candidate biomarkers for clinical validation [9]. However, most of previous studies were cross-sectional studies, and there is still insufficient evidence on the role of urinary mRNAs in predicting the progression of DN.
Weighted gene co-expression network analysis (WGCNA), a systems biology algorithm, can identify gene modules that are co-expressed and correlated with clinical traits, providing us a powerful tool to identify candidate biomarkers in multiple diseases [10,11]. Urinary single-cell RNA sequencing indicated the similarities in gene expression between cells of kidney and cells in urine [12]. It is suggested that the changes of gene expression in urinary cells can reflect the changes of renal tissue to some extent. In this study, we recognized the gene modules negatively related to eGFR based on available published microarray data of renal biopsies by WGCNA. Also, gene ontology (GO) and pathway enrichment analysis were conducted and candidate hub genes were identified by bioinformatics tools. Furthermore, the mRNA levels of candidate genes were detected in urine sediment and their correlations with albuminuria, eGFR and renal interstitial fibrosis were performed in validation cohorts. Finally, the patients were followed up to investigate the role of candidate biomarkers in predicting DN progression.

Gene expression data and pre-processing
Raw CEL files of microarray datasets GSE30528 and GSE30529 [13] were downloaded from the GEO Database (https://www. ncbi.nlm.nih.gov/geo/), containing mRNA data from glomeruli samples of 9 patients with DN and 13 healthy donors and tubule samples of 10 DN patients and 12 healthy donors, respectively. Raw data were pre-processed using R version 3.6.1 (R Foundation for Statistical Computing, Vienna, Austria) and Bioconductor (version 2.14; www.bioconductor.org). Background correction and normalization were performed by the robust multichip average (RMA) function in the affy (version 1.52.0) R package. Then we mapped the array probes to the respective Entrez gene symbol using array annotation data in the hgu133plus2.db (version 3.11) R package. Clinical data of GSE30528 and GSE30529 were downloaded from the Nephroseq platform (https://www.nephroseq.org/), including age, sex, ethnicity and eGFR.

Bioinformatics analysis
According to the protocols of the WGCNA R package (version 1.70.3), gene co-expression modules were created based on transcript and clinical data for GSE30528 and GSE30529. Module eigengenes (MEs) were calculated to summarize and represent the gene modules. Then the MEs were analyzed for correlations with clinical traits including age, sex, ethnicity and eGFR.
To explore whether module genes share a similar biological function, the visualization of GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis for the selected gene modules were performed using the R package clusterProfiler (version 3.11). Adjusted P-values <0.05 were considered as statistical significance. Search Tool for the Retrieval of Interacting Genes (STRING) is an online platform that collects comprehensive information of experimental and predictive interactions of proteins. It was applied to establish protein-protein interaction (PPI) networks and Cytoscape software (version 3.7.2) was used for visualization. The Maximal Clique Centrality (MCC) algorithm in cyto-hubba, a plug-in of Cytoscape software, was performed to identify hub genes in the PPI networks. In order to compare the levels of interested hub genes between DN and normal subjects, their mRNA levels were obtained from the microarray data in the two datasets. Pearson's correlation was performed between the expression levels of hub genes and eGFR.

Study participants
The study protocol was approved by the Ethical Committee of Zhongda Hospital Southeast University (2019ZDSYLL057-P01) and all of the subjects gave their written informed consent in compliance with the Helsinki Declaration. A total of 91 patients with renal biopsy-confirmed DN, 60 T2DM patients with no signs of renal damage [urinary albumin: creatinine ratio (uACR) <30 mg/g and eGFR ≥60 mL/min/1.73 m 2 ] and 61 healthy controls (HCs) were recruited from August 2018 to July 2019 from three hospitals in China (Zhongda Hospital Southeast University, First Affiliated Hospital of Zhengzhou University and Second Xiangya Hospital of Central South University). T2DM was diagnosed based on the 2021 criteria of the American Diabetes Association [14]. DN was diagnosed by renal biopsy according to the pathological characteristics of diagnostic criteria [15]. HCs were enrolled if they were non-diabetic and without clinical or laboratory evidence of kidney diseases. The exclusion criteria were as follows: infection or cancer, signs or symptoms of nondiabetic kidney disease, eGFR <15 mL/min/1.73 m 2 or receiving renal replacement therapy and suffering from a mental illness.
The glomerular lesions of DN were assessed according to the Tervaert classification [16], calculating the percentage of glomerulosclerosis. Ten randomly selected visual areas were assessed in each patient, and renal interstitial fibrosis was semi-quantitatively graded into mild (<25%), moderate (>25% but <50%) and severe (>50%) by calculating the average

Collection of urine samples and total RNA extraction
A whole-stream early-morning urine specimen for each subject was collected in a sterile container and then treated within 1 h after collection. An aliquot of 50 mL was centrifuged at 3000 g for 30 min at 4°C to obtain the urine sediments. The supernatant was then discarded and the urinary pellet was washed twice in 1 mL phosphate-buffered saline (PBS) and centrifuged at 12 500 g for 5 min at 4°C. Total RNA was isolated with the miRNeasy mini kit (217 004; Qiagen, Düsseldorf, Germany) according to the manufacturer's instructions and stored at −80°C until use.

Polymerase chain reaction (PCR) assay
The quantification of the candidate genes was performed using a one-step Taqman RT-qPCR MasterMix kit (Ace105, AceBioX, Shanghai, China). Beta-2 microglobulin (B2M) was used as the housekeeping gene. Sequences of primers and Taqman probes used in the PCR assay are listed in Table 1. Briefly, 100 ng of RNA was used in 10 μL reaction mixtures with gene-specific primers and probes according to the manufacturers protocol. Cycling procedures were set as follows: 42°C for 5 min, 95°C for 10 min, followed by 45 cycles of 95°C for 20 s and 60°C for 1 min. PCR assays were operated using an ABI PRISM7700 platform (Applied Biosystems, Waltham, MA, USA). The relative levels of candidate mRNA were normalized by B2M, and the relative expression for each gene was calculated by the comparative method 2 -CT .

Follow-up
After renal biopsy, 91 DN patients were followed up for at least 1.5 years until December 2020. The clinical strategy was made by the nephrologist alone and was not influenced by the study. Renal function tests, including serum creatinine and blood urea, were examined at least every 3 months. As mentioned, eGFR was calculated by the Chronic Kidney Disease Epidemiology Collaboration equation according to the concentration of serum creatinine. ESRD (eGFR <15 mL/min/1.73 m 2 or receiving renal replacement therapy) was defined as the primary endpoint. The overall design of our study is summarized in the Supplementary data, Figure S1.

Statistical analysis
Statistical analysis was performed using R version 3.6.1 SPSS Statistics 25 (IBM, Armonk, NY, USA) and GraphPad Prism 7.00 (GraphPad Software, San Diego, CA, USA). The data normality was determined by the Shapiro-Wilk test. The normally distributed variables were given as mean ± standard deviation (SD) and non-normally distributed variables were given as median with interquartile range (IQR). Depending on the distribution type of values, unpaired t-test or Mann-Whitney U-test was used to evaluate the significant differences between DN patients and HCs; one-way ANOVA analysis of variance or Kruskal-Wallis test was performed for comparisons of data with more than two groups and Holm-Sidak's or Dunn's multiple comparisons test was performed in groups with statistical significance. Pearson's or Spearman's correlation analysis was used to evaluate the correlation between gene expression levels and eGFR or uACR for each subject. To determine the independent risk factors of the primary endpoint, multivariate Cox proportional hazards regression was developed to calculate the hazard ratio (HR) and 95% confidence interval (95% CI) and variate with P-values <0.1 in the univariate analysis were added to the regression model. All P-values were two tailed and P < 0.05 was considered a statistically significant difference.

Gene modules negatively associated with eGFR identified by WGCNA mainly enriched in regulation of chemokine signaling pathway
In order to find the key modules negatively associated with eGFR, WGCNA was performed on the microarray data and clinical information of GSE30528 and GSE30529. We found 19 modules in glomeruli samples and 14 modules in tubule samples. MEs were calculated as the representative for modules, then Pearson's correlation coefficients were evaluated between MEs and the clinical traits of eGFR. According to the heatmap of module-trait correlations, we identified that the steel blue module was the highest negatively correlated with eGFR (r = −0.70, P = 0.008) in glomeruli samples ( Figure 1A), including 515 genes. To reveal potential biological functions of the genes within the steel blue modules, we conducted GO and KEGG pathway enrichment analysis. In glomeruli samples, the top 10 GO terms were significantly enriched in T-cell activation, regulation of lymphocyte activation, immune response and leukocyte migration ( Figure 1B). According to KEGG analysis, genes were significantly enriched in the chemokine signal pathway, phagosome, cell adhesion and T cell receptor signal pathway ( Figure 1C). As shown in Supplementary data, Figure S2A, the PPI network was constructed using STRING. In tubule samples, the yellowgreen module showed the highest negative correlation with eGFR (r = −0.78, P = 2e-05) (Figure 2A), including 525 genes. The top 10 GO terms were significantly enriched in extracellular matrix, regulation of cytokines secretion and regulation of immune system process ( Figure 2B). According to the KEGG enrichment analysis, genes were significantly enriched in tuberculosis, extracellular matrix receptor interaction and phagosome ( Figure 2C). As shown in Supplementary data, Figure S2B, the PPI network was constructed using STRING.
To compare the expression of chemokines screened in these two datasets, mRNA levels of CCL5, CXCL8 in GSE30528 and CXCL1, CXCL6 and CXCL12 in GSE30529 were extracted from the uploaded data. The log2-transformed mRNA level of CXCL8 was not statistically different between the two groups (Supplementary data, Figure S3). As shown in Figure 3A-D, the mRNA levels of CCL5, CXCL1, CXCL6 and CXCL12 were significantly higher in DN than in the HC group. To evaluate the potential roles of candidate chemokines in detecting the progression of DN, Pearson's correlation analysis between log2-transformed mRNA levels and eGFR were performed. As shown in Figure 3E-H, the expression levels of CCL5, CXCL1, CXCL6 and CXCL12 were negatively correlated with eGFR.

Characteristics of subjects
A total of 212 subjects were enrolled in the study (DN, n = 91; T2DM, n = 60; HCs, n = 61). The baseline clinical, laboratory and pathological characteristics of the subjects are summarized in Table 2. There was no significant difference in sex ratios among three groups. The HCs had normal albuminuria, blood glucose levels and eGFR. Patients with T2DM were older and had higher levels of glycosylated hemoglobin than the other two groups. Patients with DN had lower levels of eGFR, plasma albumin and hemoglobin, and higher levels of albuminuria and uric acid than patients with T2DM.

Urinary CCL5 and CXCL1 mRNAs were upregulated in the DN patients
The expression levels of candidate genes (CCL5, CXCL1, CXCL6 and CXCL12) were validated in the urinary sediment in cohorts by one-step Taqman PCR assay. B2M and the four chemokines could be amplified in the samples collected from all the subjects. As shown in Table 2 and Figure 4, there was no significant difference in urinary expression of CXCL6 and CXCL12 among three groups and the expression levels of CCL5 and CXCL1 were significantly upregulated in the DN group. Furthermore, CCL5 was also positively correlated with uACR (r = 0.2693, P = 0.0098) and negatively correlated with eGFR (r = −0.4388, P < 0.0001); CXCL1 was negatively correlated with eGFR (r = −0.2275, P = 0.0301) ( Figure 5). As shown in Figure 6, there was no significant difference in the expression of CCL5 and CXCL1 in different Tervaert stages of DN. In the group of severe renal interstitial fibrosis, the expression levels of urinary CCL5 (P = 0.0274) and CXCL1 (P = 0.0438) were higher than those in the mild and moderate renal interstitial fibrosis groups. There was no statistical difference for the expression levels of CXCL6 and CXCL12 in   different Tervaert stages and grades of renal interstitial fibrosis (Supplementary data, Figure S4).

Risk factors to ESRD
After a median follow-up of 21 months, 38 of 91 DN patients progressed to ESRD and 4 patients died from complications associated with ESRD. During the study period, six patients were lost to follow-up. Clinical parameters, histological damage and mRNA levels of candidate genes that had P-values <0.1 in univariate Cox regression analysis were entered into the prediction model. As shown in Table 3, univariate Cox regression analysis revealed that hemoglobin, glycosylated hemoglobin, albumin, eGFR, uACR, percentage of glomerulosclerosis, proportion of tubulointerstitial fibrosis and urinary mRNA level of CCL5 were selected to develop the prediction model. The results of the multivariate Cox proportional hazards regression for ESRD are presented in Table 4

DISCUSSION
As an algorithm for gene co-expression network construction, WGCNA can divide the whole transcriptome into several modules according to the similarity of gene expression patterns and analyze the correlation between these modules and specific clinical characteristics. Thus it might be a powerful tool for finding the gene sets closely related to eGFR and to discover the hub genes potentially related to the progression of DN. At present there have been only a few published reports applying WGCNA in DN research. In 2012, Tang et al. [17] analyzed gene expres-sion profiles and identified several transcription factors of DN by WGCNA. Shang et al. [18] identified three key long noncoding RNAs (NR_130 134.1, NR_038 335.1 and NR_029 395.1) and found they might play important roles in DN by regulating the biological process of inflammation and metabolism. More recently, a study using WGCNA suggested that HIF1A is involved in renal interstitial fibrosis in DN [19]. To our knowledge, biomarkers correlated with a decline in eGFR of DN identified by WGCNA have not yet been reported previously. The pathogenesis of DN remains complex and involves numerous distinct pathways. Previous studies have proved urinary mRNAs of podocyte markers, such as nephrin and podocin, are correlated with a decline in renal function in DN [20,21]. Furthermore, measuring the mRNA levels of podocyte markers in urine may be valuable for monitoring treatment response [22]. In the present study we analyzed the gene clusters negatively related to eGFR in microdissected glomeruli and tubules of patients with DN by WGCNA, avoiding the bias in selecting candidate biomarkers. On the one hand, through functional annotation and pathway enrichment analysis of gene modules, module genes were mainly involved in chemokine signal and inflammatory cell activation. Most of the top 10 MCC score hub genes in glomeruli and tubule samples were members of the chemokines or chemokine receptors. It suggested that chemokines may play an important role in the progression of DN. On the other hand, the candidate chemokines were identified by hub gene screening and differential expression analysis. After multicenter clinical verification by quantitative PCR, it was found that urinary CCL5 and CXCL1 were related to the decline in renal function and the degree of renal interstitial fibrosis, and CCL5 could independently predict the disease progression of DN. Our study indicated that chemokines play an important role in the progression of DN and identified urinary CCL5 mRNA as a novel prognostic biomarker for DN, providing a novel approach of systems biology to study DN biomarkers. Renal interstitial fibrosis was semi-quantitatively graded into mild (<25%), moderate (>25% but <50%) and severe (>50%). RIF: renal interstitial fibrosis. P < 0.05 was considered statistically significant. Traditionally, DN has long been regarded as a noninflammatory glomerular lesion induced by hemodynamic and metabolic disorders. Increasing evidence from experimental and clinical studies has demonstrated that both systemic and local inflammation play crucial roles in the occurrence and progression of DN [23]. Consequently, inflammatory cytokines have attracted investigators' attention as potential biomarkers of DN. Elevated concentrations of inflammatory cytokines, including tumor necrosis factor (TNF)-α and TNF receptors (TNFRs), interleukin (IL)-1, IL-18 and IL-6 have been detected in the peripheral blood of patients with DN, with some increasing with disease progression [24]. For example, numerous prospective studies have suggested a strong association between plasma TNF-R1/2 and increased risk of DN progression, including a decline in renal function, occurrence of ESRD and all-cause mortality [25][26][27]. Chemokines, together with their receptors, induce inflammation by promoting activation and recruitment of circulating leukocytes toward the kidney. Macrophages are considered as key inflammatory cells involved in kidney repair or fibrosis [28]. In DN patients, the interstitial infiltration of inflammatory cells, such as macrophages and lymphocytes, could be attracted by chemokines released from injured renal tissue [29]. CCL2 is one of the most widely studied chemokines in renal injury [30]. The evidence of urinary CCL2 as a potential biomarker of detecting progression [31][32][33] and an intervention target [34] of DN demonstrated the potential application of chemokines in clinical utility. A proteomics study of patients with T1DM or T2DM identified a panel of inflammatory signatures of kidney outcome risk, including 17 circulating inflammatory cytokines (including CCL14 and CCL15) that were associated to progression to ESRD followed for 8-13 years [35]. It also emphasized the critical roles of chemokines in the progression of DN.
In the present study, we found that the levels of CXCL1 and CCL5 in urine were increased in DN, which may be associated with the deterioration of renal function and renal interstitial fibrosis. Several reports have been published regarding the role of these two chemokines in experimental and clinical studies of DN. A previous study indicated CXCL1 was increased in glomerular podocytes of OVE, db/db and STZ-induced DM mice [36]. However, there is not sufficient evidence for the expression of CXCL1 in patients with DN. CCL5 is a potent chemoattractant for macrophages, T cells and granulocytes and is expressed by multiple cell types, including mesangial cells, fibroblasts and tubular epithelial cells in the kidney. In STZinduced DM mice, CCL5 was increased with inflammation and fibrosis during kidney injury [37]. Moreover, upregulation of CCL5 was observed in renal biopsy tissues from patients with DN, mainly located in tubular cells, and expression of CCL5 has a positive correlation with proteinuria and interstitial inflammation [38]. Increased plasma or urinary levels of CCL5 have been reported in scattered studies in patients with T2DM and were associated with the development of microalbuminuria [39,40]. The present study found that urinary CCL5 and CXCL1 mRNA levels could reflect the decline of renal function in patients with DN. In addition, urinary CCL5 mRNA could potentially predict the progression to ESRD. These findings would provide novel methods for assessing the prognosis of DN and warrant further clinical investigation. Some shortcomings exist in this study, such as the sample size was relatively small and the follow-up period was also short. Clearly, further studies with larger samples and a longer follow-up period to confirm the role of CCL5 mRNA in the progression of DN are needed.

CONCLUSIONS
To sum up, through WGCNA analysis, we found that numerous chemokine members might be related to the progression of DN. Through clinical verification and prospective follow-up, it was found that the expression of urinary CCL5 mRNA and CXCL1 mRNA is associated with lower eGFR and severe renal interstitial fibrosis, and urinary CCL5 mRNA could be an independent risk factor for DN progression to ESRD. This study may provide a methodological reference for identifying key molecules related to specific phenotypes in DN.

DATA AVAILABILITY STATEMENT
The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.