Gene expression profiles of GSE73680 were downloaded from the Gene Expression Omnibus (GEO) database, and expression profiling arrays were generated using GPL17077. Biopsies were taken from human renal papillary tip tissues during endoscopic kidney stone surgery. In total, 24 samples were extracted in our study, including 13 CaOx plaque tissues, 5 CaP plaque tissues, and 6 normal control tissues. Then, the following steps were sequentially applied: removing the control probe, removing the probes that matched no gene symbol, removing the gene symbol with LOC or MT, and averaging duplicated gene probes. Finally, 28,892 genes were used for further analysis.

To characterize DEGs, the limma package in the R software (version 4.2.2) was used to identify DEGs between the CaOx group and the control group with the threshold of p-value < 0.05 and |log2 fold change (FC)| >1. The limma package was also used to identify DEGs between the CaP group and the control group with the threshold of p-value < 0.05 and | FC| >1. The Venn diagram was constructed using the Venn diagram package. The volcano plot of DEGs was drawn by the ggplot2 package, and the heatmap of DEGs was generated using the pheatmap package in the R software.

To identify the potential functions and pathways between these DEGs, Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted using David, and results were visualized using R ggplot2 package. The GO terms of biological processes (BP), molecular functions (MF), and cellular components (CC) were respectively evaluated. The top ten terms were visualized if there were more than ten terms.

To systematically analyze the biological functions of the obtained DEGs, the DEGs identified previously were mapped into the online search tool STRING database, which predicted the protein functional associations and PPI. Then, the Cytoscape software (version 3.7.2) was used for constructing and visualizing the transcriptional regulatory network of common DEGs. To identify hub genes, the sides of each node were estimated and the genes were sorted based on the rank of the connection number of each gene.

In our study, CCL7 was finally identified as the top hub gene for GSEA. Thirteen CaOx plaque tissues and 5 CaP plaque tissues were divided into high- and low-expression groups by the median expression value of CCL7. Then, GSEA analysis was applied to compare the differences in enhanced functions or pathways between high- and low-CCL7 expression groups using clusterprofiler package in the R software. Adjusted p-value < 0.05 was considered statistically significant.

In total, 9 patients with primarily CaOx mixed with CaP stones and 9 healthy subjects with similar gender and age were enrolled from April to October 2023. Patients’ kidney stone samples were collected after percutaneous nephrolithotomy, and stone composition was confirmed by the infrared spectrum analysis method. Patients with any evidence of kidney infection or dysfunction, a history of urinary tumors, or those younger than 18 years were excluded. On admission, midstream morning urine was collected from all participants before giving any medicine. This study was approved by the Human Ethics Committee of the Lishui Central Hospital (approved number 2023054). All patients signed the informed consent before participation.

Fresh urine was centrifuged 5000 rpm for 5 min, then discarded the supernatant. Total RNA was extracted from urine using RNAUrineKit (Solarbio, R1300-100T, China) according to the manufacturer’s instruction. The quality of all RNA samples was checked before downstream molecular biology applications. We used a tool for RNA quality assessment, called the RNA Integrity Number (RIN) which developed by Agilent Technologies. The RIN ranging from1 (degraded) to 10 (intact) was established by Agilent Bioanalyzer 2,100(Agilent Technologies). A RIN higher than five was used for downstream application. The concentration was assessed by NanoDrop 2000 (Thermo Fisher, United States). Extracted RNA samples with more than 50 ng/μL concentrations were used for analysis.

The reverse transcription kit (Takara, RR037A, Japan) was applied to synthesize cDNA (CCL7, OR10A5, OR9A2, and OR1L3) according to the manufacturer’s protocol. RT-PCR was accomplished using HotStart™ 2X SYBR Green qPCR Master Mix (APExBIO, K1070, USA) to amplify these genes with primer pairs (Table 1). Total RNA of 1 μg was used for cDNA synthesis. All primers were synthesized by Sangon Biotech. formula 2^-△△Ct was used to analyze the relative changes in mRNA expression from RT-PCR (Pfaffl, 2001).

Statistical analysis was performed through R software (version 4.2.2) and GraphPad Prism software (version 9, San Diego, CA). Continuous normally distributed variables are presented in mean ± standard deviation (SD), and experimental results are expressed as mean ± standard error (SEM). Student’s t-test was applied to compare the stone and control group samples, with p < 0.05 showing a statistically significant difference. Receiver Operating Characteristic Curve (ROC) was constructed and used to evaluate the diagnostic value of CCL7 for patients with KSD.

In total, 2,885 and 506 DEGs were identified from CaOx vs. control samples and CaP vs. control samples, respectively (see Supplementary Tables S1, S2). Volcano plots showed the distribution of these DEGs (Figures 1A, B). Subsequently, the intersection of these DEGs was identified using the Venn diagram, containing 155 significantly co-up-regulated genes and 64 co-down-regulated genes (Figure 1C and Supplementary Table S3). Therefore, the number of co-up-regulated DEGs was more than the number of co-down-regulated DEGs. The 30 most upregulated or downregulated DEGs were visualized in a heatmap (Figures 2A, B).

To investigate the functions and mechanisms of these DEGs, GO and KEGG pathway enrichment analyses of upregulated and downregulated genes were performed. According to the results of enrichment analysis, 18 GO terms of DEGs (p < 0.05), including 6 BP, 6 CC, and 6 MF, were displayed (Figure 3A). The GO analysis showed that DEGs were significantly enriched in detection of chemical stimulus in sensory perception, detection of chemical stimulus in sensory perception of smell, and olfactory receptor activity (see Supplementary Table S4). Moreover, 6 KEGG pathways were over-represented in the DEGs, and the olfactory transduction pathway was significantly enriched. The results of the KEGG enrichment analysis are shown in Figure 3B and Supplementary Table S5.

To investigate the relationships among these DEGs at the protein level, the PPI of the DEGs was obtained using STRING with a confidence score>0.4. Then, the PPI network of these DEGs was visualized using Cytoscape. Figure 4A presented a PPI network with 13 nodes and 27 edges, including 10 upregulated genes and 3 downregulated genes. In total, 10 genes were identified as hub genes according to the cytoHubba plugin of Cytoscape. They were all upregulated genes, and CCL7 was the top hub gene. The gene symbols and the detailed nodes and edges of the PPI network are shown in Figure 4B.

Eighteen RP papillary tissues were divided into two subgroups based on the median expression of CCL7. The olfactory transduction, maturity-onset diabetes of the young, linoleic acid metabolism, and fat digestion and absorption were significantly enriched in the high-CCL7 subgroup (Figure 5).

There was no significant difference in sex (both M/F 5:4) and age (52.8 ± 10.5 vs 54.2 ± 11.4 years) between the stone group and the control group (Table 1). The results of RT-PCR showed that CCL7 level in the stone group was significantly higher than that in the control group, which was consistent with the results of bioinformatic analysis. For the olfactory transduction pathway, the expression of OR10A5, OR9A2, and OR1L3 was significantly upregulated in the stone group compared with the control group (Figure 6), which was also consistent with bioinformatic results.

The value of CCL7 in both renal papillary tip tissues (24 patients) and urine samples from our PCR results (18 patients) in diagnosing KSD was investigated by plotting a ROC curve. As shown in Figures 7A, B, CCL7 in renal papillary tip tissues predicted the extent of KSD with AUC = 0.935 (p < 0.001), and CCL7 in urine samples predicted the extent of KSD with AUC = 0.985 (p < 0.001). They both had high sensitivity and a specificity.