- 1Department of Urology, Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
- 2Division of Cardiology, Department of Internal Medicine, Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
Background: The pathogenesis of urolithiasis remains unclear, making the development of medications for treatment and prevention stagnant. Randall’s plaques (RPs) begin as interstitial calcium phosphate crystal deposits, grow outward and breach the renal papillary surface, acting as attachment for CaOx stones. Since matrix metalloproteinases (MMPs) can degrade all components of extracellular matrix (ECM), they might participate in the breach of RPs. Besides, MMPs can modulate the immune response and inflammation, which were confirmed to be involved in urolithiasis. We aimed to investigate the role of MMPs in the development of RPs and stone formation.
Methods: The public dataset GSE73680 was mined to identify differentially expressed MMPs (DEMMPs) between normal tissues and RPs. WGCNA and three machine learning algorithms were performed to screen the hub DEMMPs. In vitro experiments were conducted for validation. Afterwards, RPs samples were classified into clusters based on the hub DEMMPs expression. Differentially expressed genes (DEGs) between clusters were identified and functional enrichment analysis and GSEA were applied to explore the biological role of DEGs. Moreover, the immune infiltration levels between clusters were evaluated by CIBERSORT and ssGSEA.
Results: Five DEMMPs, including MMP1, MMP3, MMP9, MMP10, and MMP12, were identified between normal tissues and RPs, and all of them were elevated in RPs. Based on WGCNA and three machine learning algorithms, all of five DEMMPs were regarded as hub DEMMPs. In vitro validation found the expression of hub DEMMPs also increased in renal tubular epithelial cells under lithogenic environment. RPs samples were divided into two clusters and cluster A exhibited higher expression of hub DEMMPs compared to cluster B. Functional enrichment analysis and GSEA found DEGs were enriched in immune-related functions and pathways. Moreover, increased infiltration of M1 macrophages and enhanced levels of inflammation were observed in cluster A by immune infiltration analysis.
Conclusion: We assumed that MMPs might participate in RPs and stone formation through ECM degradation and macrophages-mediated immune response and inflammation. Our findings offer a novel perspective on the role of MMPs in immunity and urolithiasis for the first time, and provide potential biomarkers to develop targets for treatment and prevention.
Introduction
Urolithiasis is a common urological disease with an increasing incidence rate and a high recurrence rate (1). Recurrent episodes of stones lead to renal injury and chronic kidney disease, posing a huge economic burden on individuals and society (2). Although rapid advances of surgical techniques have greatly improved the stone removal efficiency, the development of medical drugs to treat and prevent urolithiasis has stalled, as its pathogenesis is not fully understood. Calcium oxalate (CaOx) stones are the most common stone type (3), and hypercalciuria and hyperoxaluria are major risk factors for stone formation (4). Randall’s plaques (RPs) theory is the most widely accepted theory about CaOx stone formation. RPs begin as calcium phosphate (CaP) crystal deposits in the renal interstitium, grow and extend outwards reaching the renal papillary surface, and become exposed to the pelvis urine. CaOx crystals could attach to the exposed sites and grow into CaOx stones (5, 6). Apparently, RPs are subepithelial deposits that must have their surfaces breached. Some researchers believe that papillary surface epithelium can be breached through the involvement of matrix metalloproteinases (MMPs) and/or the sheer force of the growing plaque (7). Still, how such a breach happens remain obscure. MMPs are a family of zinc and calcium dependent endopeptidases which are able to degrade all components of the extracellular matrix (ECM) and process various bioactive substances (8). Currently, the MMP family consists of 23 members in human, which can be divided into 6 major types, including collagenases, gelatinases, stromelysins, matrilysins, membrane-type MMPs, and other MMPs (8). Under physiological conditions, MMPs participate in tissue repair and remodelling, morphogenesis, wound healing, cell proliferation and migration (9, 10). Dysregulation of MMPs might result in tissue destruction and homeostasis disorders and contribute to various diseases, such as cancer, atherosclerosis, aortic aneurysms, arthritis, tissue ulcers and fibrosis (9, 10). In addition, MMPs have been implicated in a wide range of kidney diseases, such as kidney fibrosis, glomerular disease, and diabetic kidney disease (11). So far, few studies have been reported the relationship between MMPs and urolithiasis.
Besides, immunity and inflammation were tightly associated with urolithiasis (6). Immune cells and pro-inflammatory molecules were shown to be increased in RP tissues and macrophages seem to participate critically in stone formation (12). M1 macrophages can induce inflammatory response to accelerate renal injury and CaOx crystals depositions, while M2 macrophages can phagocytose and degrade crystals (13). MMPs were also found to be involved in immunity and inflammation, serving as modulators of immune response and inflammatory processes, which can influence and shaped the phenotype of immune cells (14). Therefore, we assumed that MMPs might play a pivotal role in the development of RPs. However, the precise mechanism is still a mystery.
In this study, we first used the dataset GSE73680 to identify differentially expressed MMPs (DEMMPs) between normal renal papillary tissues and RPs. Then, we performed weighted gene co-expression network analysis (WGCNA) and three machine learning algorithms to screen the hub DEMMPs. Afterwards, we divided RPs samples into two clusters based on the hub DEMMPs expression pattern, identified differentially expressed genes (DEGs) between them and conducted further functional enrichment analysis and gene set enrichment analysis (GSEA). Moreover, we evaluate the immune infiltration levels between two clusters. We believe our findings will offer a novel insight into the role of MMPs in urolithiasis and provide potential biomarkers to develop therapeutic targets for the disease.
Materials and methods
Acquisition and processing of data
The flowchart of the study is shown in Figure 1. The dataset GSE73680 was obtained from gene expression omnibus database (GEO, http://www.ncbi.nlm.nih.gov/geo/), which contains 6 normal renal papillary tissues from non-stone formers and 29 RPs from stone formers (12). The annotation information provided in the platform GPL17077 were used to transform the probe IDs into gene symbols.
Identification of DEMMPs
The MMP genes were obtained through HUGO Gene Nomenclature Committee (https://www.genenames.org/). The Wilcox test was applied to identify DEMMPs between normal tissues and RPs, and the results were shown as a boxplot using the “ggpubr” package. Besides, the correlations between MMPs were calculated by Spearman correlation analysis and the correlation matrices were visualized using the “corrplot” package. The protein-protein interaction (PPI) network of MMPs was established by the STRING database (https://string-db.org/) with minimum interaction score ≥ 0.4. We then calculate the sides of each node and sort genes based on the rank of the connection number of each gene.
Weighted gene co-expression network analysis
The “WGCNA” package was utilized to construct the co-expression networks for identification of hub module related with urolithiasis (15). Initially, the Pearson’s correlation matrices were performed for all pair-wise genes and a weighted adjacency matrix was constructed with the formula amn = |cmn|β (cmn = Pearson’s correlation between gene m and gene n; amn = adjacency between gene m and gene n). Then, the soft threshold power value, parameter β was calculated, which emphasized strong correlations between genes and penalize weak correlations. A suitable parameter β was applied to construct a scale-free network. Next, the weighted adjacency matrix was transformed into a topological overlap measure (TOM) matrix, which could measure the network connectivity of a gene defined as the sum of its adjacency with all other genes for the network generation (16). Afterwards, average linkage hierarchical clustering was performed to classify genes with similar expression profiles into the same gene modules based on the TOM‐based dissimilarity measure with a minimum gene group size of 50 for the genes dendrogram (17). Finally, we calculated the gene significance (GS), module membership (MM), and associated modules with clinical traits to identify the urolithiasis-related modules with the highest coefficient square (R2) and the P-value < 0.05. We took the intersection between DEMMPs and genes in urolithiasis-related modules to obtain the hub DEMMPs
Machine learning
Three machine learning algorithms were used to screen the hub DEMMPs. The least absolute shrinkage and selection operator (LASSO) is a regression approach used for dimensionality reduction and feature selection to increase predicted accuracy and model comprehensibility (18), which was performed by the “glmnet” package (19). Support vector machine recursive feature elimination (SVM-RFE) is a powerful method to rank features and to select the significant ones for classification (20), which was performed by the “kernlab” and “e1071” packages (21). Random forest (RF) is an ensemble learning method for classification by constructing a multitude of decision trees at training time (22), which was performed by the “randomForest” package (23). DEMMPs identified more than three times in the above four analysis (WGCNA, LASSO, SVM-RFE, RF) were considered as hub DEMMPs. The receiver operating characteristic (ROC) was constructed and the area under the ROC curve (AUC) was calculated to evaluate the capacity of hub DEMMPs to discriminate RPs from normal tissues using the “pROC” package (24).
Cell experiments and validation of hub DEMMPs
The normal rat kidney proximal tubular epithelial cell line (NRK-52E) was obtained from the Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China). Cells were cultured in Dulbecco’s modified eagle medium (DMEM; Hyclone, USA) with 10% fetal bovine serum (FBS; Gibco, USA) at 37°C in a humidified atmosphere with 5% CO2. The normal human kidney proximal tubular epithelial cell line (HK-2) was purchased from Procell (Wuhan, China) and cultured in DMEM-F12 medium (Hyclone, USA) with 10% FBS at 37°C in a humidified atmosphere with 5% CO2. NRK-52E and HK-2 cells were treated with oxalate (1 mmol/L) for 24 h, and cells without treatments were used as a control group.
Total RNA was extracted using a TRIzol reagent (Vazyme, Nanjing, China). Then, the purified RNA was reverse-transcribed into cDNA using a cDNA synthesis kit (Yeasen, Shanghai, China). Next, real-time quantitative polymerase chain reaction (RT-qPCR) was performed using SYBR Green Master Mix reagent (Yeasen, Shanghai, China) on the PCR System. The mRNA expression levels of genes were calculated using the 2-ΔΔCt method. The sequences of primers used in this study were presented in Table 1.
Consensus clustering analysis
We used the “ConsensusClusterPlus” R package (25) to conduct an unsupervised hierarchical clustering analysis on 29 RPs samples based on the mRNA expression of hub DEMMPs. We used consensus matrix plots, consensus cumulative distribution function (CDF) plots, and the relative change in area under the CDF curve to determine the optimal number of clusters. After that, principal component analysis (PCA) between clusters was performed to observe the differences between clusters. The “limma” package was then used to identify DEGs between clusters with the threshold of |log2(fold change)|>2 and the P-value < 0.05. DEGs were visualized as a volcano plot and a heatmap.
Functional enrichment analysis and GSEA
To further predict the biological functions of DEGs, Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using the “clusterProfiler” package (26). The adjusted P-value < 0.05 was the threshold. We also applied GSEA to explore the significant functional differences between clusters using the “clusterProfiler” package. The gene sets “c5.go.v2022.1.Hs.symbols.gmt” and “c2.cp.kegg.v2022.1.Hs.symbols.gmt” were obtained from the GSEA website MsigDB database (http://www.gsea-msigdb.org/gsea/msigdb/). Gene sets were considered as significantly enriched by False discovery rate (FDR) < 0.25, the adjusted P-value < 0.05 and |normalized enrichment score (NES)| > 1.
Immune infiltration analysis
We used two methods to investigate the immune infiltration levels between two clusters. CIBERSORT (https://cibersortx.stanford.edu/) was applied to evaluate infiltration status of 22 distinct immune cells in two clusters, and the outcomes were visualized as a violin plot using the “vioplot” package. The correlations between the hub DEMMPs and immune cells were evaluated by Spearman correlation analysis and were visualized as lollipop plots. Single sample gene set enrichment analysis (ssGSEA) was conducted to quantify immune infiltration levels by calculating the enrichment scores of 16 immune cells and 13 immune-related functions in two clusters using the “gsva” package (27, 28), and the results were visualized as a boxplot using the “ggpubr” package.
Statistical analysis
Data analysis and statistical analysis were performed by R software (version 4.1.3) and GraphPad Prism (version 8.0). The statistical difference between two groups were assessed by the Wilcox test or the unpaired t-test. The correlations between the genes were analyzed by a spearman correlation test. The P-value < 0.05 was considered statistically significant.
Results
Identification of DEMMPs
We first screened DEMMPs between normal tissues and RPs using the Wilcox test. A total of five MMPs, including MMP1, MMP3, MMP9, MMP10, MMP12 were identified as DEMMPs, and all of them were significantly upregulated in RPs (Figures 2A, B). Through the correlation analysis, we found all DEMMPs were positively correlated with each other (Figure 2C). The PPI network was conducted to investigate the interactions between MMPs at the protein level (Figures 2D, E). For DEMMPs, MMP9 have the highest numbers of interacted proteins, followed by MMP1, MMP3, MMP10, and MMP12. Besides, the interactions between DEMMPs were shown in Figure 2F.
Figure 2 Identification of DEMMPs. (A) The boxplot of 26 MMPs expression in normal tissues and RPs. (B) The boxplot of five DEMMPs in normal tissues and RPs. (C) Correlation analysis between 26 MMPs. (D) PPI networks of 26 MMPs. (E) The histogram showing the number of adjacent nodes. (F) PPI network of five DEMMPs. * represents P < 0.05, ** represents P < 0.01; *** represents P < 0.001.
Screening of hub DEMMPs via WGCNA
We performed WGCNA to explore modules associated with urolithiasis. The samples were clustered hierarchically to remove outliers, and no outliers were detected and removed in this study (Figure 3A). Next, we chose β = 8 as the soft threshold power value to make the scale-free R2 reach 0.8 and construct a scale-free network (Figure 3B). A hierarchical clustering tree was constructed and modules were detected by dynamic tree cutting, and modules with highly correlated eigengenes were merged (minimal module size = 50, merge height = 0.25) (Figure 3C). A total of 14 color modules were determined (Figure 3D). Among them, the royalblue module (r = 0.41, P = 0.01), the green module (r = 0.39, P = 0.02) and the tan module (r = 0.34, P = 0.04) showed significantly positive correlations with urolithiasis (Figure 3E). Subsequently, we extracted DMMMPs and genes in urolithiasis-related modules for intersection and found MMP1, MMP3, MMP9, MMP12 were in the royalblue module and MMP10 were in the green module, suggesting all of DEMMPs might be implicated in RPs and urolithiasis (Figure 3F).
Figure 3 Screening of hub DEMMPs via WGCNA. (A) Hierarchical clustering tree of 6 normal tissues and 29 RPs gene expression patterns. (B) Identification of power value. The red line represents R2 > 0.8 when the power value β is 8. (C) Module eigengene dendrogram presented the relationship of the modules generated by the clustering analysis. (D) Clustering dendrogram and merging of the gene co-expression modules. Each color represents one module. (E) The heatmap of the correlation between modules and clinical traits. The correlation coefficient and P value between the module and clinical traits are shown at the row-column intersection. (F) The Venn plot showing the intersection between DEMMPs and the genes in urolithiasis-related modules.
Screening of hub DEMMPs via machine learning
We further performed three machine learning algorithms using the expression profiles 26 MMPs to identify potential MMPs associated with RPs and urolithiasis. For LASSO regression analysis, six variables, MMP1, MMP3, MMP9, MMP10, MMP27 and, MMP28, were screened as diagnostic markers for urolithiasis (Figure 4A). For the SVM-RFE algorithm, the error was minimized when the number of variables was six, including MMP10, MMP17, MMP9, MMP12, MMP1, and MMP3 (Figure 4B). For RF algorithm, six candidate genes with relative importance scores greater than 0.5 were identified, including MMP10, MMP9, MMP12, MMP17, MMP1, and MMP3 (Figure 4C). Combined with the outcomes of WGCNA, we found all of DEMMPs were identified to be associated with urolithiasis by at least three out of four approaches used. Therefore, we considered all of DEMMPs as hub DEMMPs for further analysis.
Figure 4 Screening of hub DEMMPs via machine learning. (A) LASSO analysis. Vertical dashed lines are plotted at the best lambda. (B) SVM-RFE algorithm for feature gene selection. (C) Ranking of the relative importance of genes in the RF classifier. (D) ROC curve of hub DEMMPs in urolithiasis diagnosis.
Next, we evaluated the diagnostic performance of hub DEMMPs by constructing ROC curves. Notably, MMP10 had the highest AUC among the seven hub genes, with a value of 0.948. The AUC values for other genes were 0.845 for MMP1, 0.828 for MMP9, 0.805 for MMP3, and 0.805 for MMP12, respectively (Figure 4D). These results indicated that all five gene signatures had powerful diagnostic values.
Validation of hub DEMMPs
We built a cell model of urolithiasis to validate hub genes. NRK-52E cells exposed to oxalate showed higher expression of MMP1 (10.30 fold, P = 0.0139), MMP3 (11.36 fold, P = 0.0008), MMP9 (4.35 fold, P = 0.0015), MMP10 (19.15 fold, P < 0.0001), and MMP12 (6.80 fold, P = 0.0002) compare to control (Figure 5A). HK-2 cells exposed to oxalate also showed higher expression of MMP1 (3.42 fold, P = 0.0054), MMP3 (13.53 fold, P < 0.0001), MMP9 (1.64 fold, P = 0.0009), MMP10 (3.81 fold, P = 0.0194), and MMP12 (2.21 fold, P = 0.0059) compare to control (Figure 5B). The results were consistent with bioinformatics screening.
Figure 5 Validation of hub DEMMPs. (A) The expression of hub DEMMPs in NRK-52E cells after exposure to oxalate for 24h. (B) The expression of hub DEMMPs in HK-2 cells after exposure to oxalate for 24h. * represents P < 0.05, ** represents P < 0.01; *** represents P < 0.001.
Identification of two clusters based on hub DEMMPs
RPs samples were clustered by the consensus clustering analysis based on the expression profiles of hub DEMMPs to discover subtypes of urolithiasis. Clustering results are relatively stable when the number of clusters was set to 2 based on the CDF curve and the CDF Delta area curve (Figures 6A, B). The consensus matrix plot showed that 29 RPs samples could be divided into two distinct clusters, i.e., cluster A (n = 11) and cluster B (n = 18) (Figure 6C). A clear distribution between cluster A and cluster B was presented in the PCA plot (Figure 6D). The expression profiles of the hub DEMMPs in the two clusters were visualized as the boxplot (Figure 6E). Cluster A showed higher expression of hub DEMMPs compared to Cluster B. Hence, cluster A was considered as high expression group and cluster B was considered as low expression group. A total of 233 DEGs were identified between cluster A and cluster B. Compared to cluster B, 204 DEGs were upregulated and 29 DEGs were downregulated in cluster A (Figures 6F, G).
Figure 6 Identification of two clusters based on hub DEMMPs. (A) The consensus cumulative distribution function (CDF) plot. (B) The relative change in area under the CDF curve. (C) The consensus matrix plot when consumption k = 2. (D) The PCA plot showing the distribution of two clusters. (E) The boxplot of five DEMMPs in cluster A and cluster B (F) The heatmap of the top 20 DEGs between two clusters. (G) The volcano plot of all DEG between two clusters. * represents P < 0.05, ** represents P < 0.01; *** represents P < 0.001.
Functional distinctions between two clusters
We further applied GO enrichment analysis and KEGG pathway analysis to explore the potential role of DEGs. Among the GO analysis results (Figure 7A), leukocyte mediated immunity, leukocyte chemotaxis, positive regulation of cell-cell adhesion, positive regulation of cytokine production, and myeloid leukocyte migration were the five most significant biological processes (BP); MHC class II protein complex, MHC protein complex, and clathrin-coated endocytic vesicle membrane were the three most significant cellular components (CC); immune receptor activity, CXCR chemokine receptor binding, and MHC class II protein complex binding were the three most significant molecular functions (MF). According to the KEGG analysis, DEGs were enriched in rheumatoid arthritis, viral protein interaction with cytokine and cytokine receptor, cytokine-cytokine receptor interaction, and so on (Figure 7B).
Figure 7 Functional enrichment analysis and GSEA between two clusters. (A) GO analysis of DEGs. (B) KEGG analysis of DEGs. (C) GSEA showing the enriched GO terms in cluster A (D) GSEA showing the enriched GO terms in cluster B (E) GSEA showing the enriched KEGG pathways in cluster A (F) GSEA showing the enriched KEGG pathways in cluster B.
We also performed GSEA to explore the enriched GO terms and KEGG pathways in two clusters. For GO terms, activation of immune response, acute inflammatory response, acute phase response, and adaptive immune response were enriched in cluster A, while carbon dioxide transport, one carbon compound transport, sensory perception of smell, olfactory receptor activity, and oxygen carrier activity were enriched in cluster B (Figures 7C, D). For KEGG pathways, allograft rejection, antigen processing and presentation, cell adhesion molecules cams, chemokine signaling pathway, and cytokine-cytokine receptor interaction were enriched in cluster A, while linoleic acid metabolism and olfactory transduction were enriched in cluster B (Figures 7E, F). Since cluster A had higher expressions of hub DEMMPs, our results suggested MMPs might be associated with immune response and inflammation.
Differences of immune infiltration between two clusters
Based on differences of functional analysis, we further investigate the relationship between immune infiltration levels between two clusters. The outcomes of CIBERSORT showed that the level of macrophages M1 was relatively high, while the level of neutrophils was relatively low in cluster A compared to cluster B (Figure 8A). The correlations between hub DEMMPs and immune cells were shown in Figure 8B. Hub DEMMPs exhibited significant positive and negative correlations with different immune cells. MMP1, MMP3, and MMP12 were positively correlated with M1 macrophages (R = 0.46, P = 0.011; R = 0.38, P = 0.044; R = 0.42, P = 0.025). MMP9 and MMP10 showed positive correlation tendency with M1 macrophages, albeit without reaching statistical significance (R = 0.21, P = 0.269; R = 0.25, P = 0.191).
Figure 8 Immune infiltration analysis between two clusters. (A) The violin plot showing different proportions of 22 immune cells between two clusters by CIBERSORT. (B) The lollipop plots showing the correlations between hub DEMMPs and immune cells. (C) The boxplot of 11 M1 macrophages markers in cluster A and cluster B (D) The boxplot showing different proportions of 16 immune cells and 13 immune-related functions between two cluster by ssGSEA. * represents P < 0.05, ** represents P < 0.01; *** represents P < 0.001.
Since cluster A showed increased M1 macrophages levels, we compared markers of M1 macrophages between 2 clusters, including TNF, IL1A, IL1B, IL6, IL12A, IL12B, NOS2, TLR2, TLR4, CD80, and CD86 (Figure 8C). The expression of NOS2, TLR2, and TLR4 were higher in cluster A compared to cluster B. All other markers tended to increase in cluster A, although not statistically significant. In addition, we compared ssGSEA scores between normal tissues and RPs and found cluster A generally had higher immune infiltration (Figure 8D). For example, for immune cells, macrophages and dendritic cells were more abundant in cluster A; for immune-related functions, inflammation-promoting and type I IFN response were enriched in cluster A.
Discussion
Development of medications for the treatment and prevention of urolithiasis has been lagging, since its pathogenesis is still unclear. It is well accepted that RPs are the origin of CaOx stones (29). RPs are initially renal interstitial CaP deposits which grow outward and across the renal papillary surface, and become exposed to pelvis urine. CaOx stones will form and develop with the attachment to RPs (5, 6). It is important to break through the renal epithelium during the development of RPs, but the exact mechanism remain an enigma. As MMPs are capable of degrading all kinds of ECM proteins, they might participate in the process of the breach to promote RPs progression. In addition, MMPs are able to modulate the immune response and inflammation, which were confirmed to be involved in urolithiasis (6). Therefore, we aimed to explore the potential role of MMPs in urolithiasis through data mining based on the public gene expression profile of RPs (GSE73680).
We first identified five DEMMPs between normal tissues and RPs, including MMP1, MMP3, MMP9, MMP10, and MMP12, and the expression of all DEMMPs was elevated in RPs. MMP1 belongs to collagenases, MMP3 and MMP10 belong to stromelysins, MMP9 belongs to gelatinases, and MMP12 belongs to other MMPs (8). Correlation analysis showed the expression of all DEMMPs were positively correlated with each other and PPI showed some interactions between them. Studies have demonstrated that MMP1 could cleaves pro-MMP9 into its active form; MMP3 can active pro-MMP1 and pro-MMP9 (10). This might be a good explanation for the results of correlation analysis and PPI.
Four approaches, including WGCNA, LASSO, SVM-RFE, and RF, were performed to screen the hub DEMMPs, and all DEMMPs were identified by at least three out of four approaches used. Therefore, we referred to all DEMMPs as hub DEMMPs. Corresponding AUCs for hub DEMMPs > 0.8 indicated they are powerful biomarkers to distinguish RPs from normal tissues. In addition, we also found that NRK-52E cells upregulated the expression of hub DEMMPs after exposure to oxalate. We assumed that MMPs expression might increase in renal tubular epithelial cells under lithogenic environment, and they would degrade ECM components to facilitate RPs to breach the renal epithelium.
In cardiovascular diseases, increased MMP1 expression in the heart led to a disruption of structural collagen and cardiac dysfunction (30). Activation of MMP3 promoted ECM damage and increased the risk for developing abdominal aortic aneurysms (31). MMP9 and MMP12 promoted vascular smooth muscle cell proliferation via beta-catenin pathway and contributed to atherosclerosis (32). MMP10 played a pivotal role in aortic valve calcification by inducing cell mineralization (33). Studies also found the levels of MMP1, MMP3, and MMP12 were positive associated with the degree of atherosclerosis and plaque stability, acting as valuable biomarkers for clinical diagnosis and prognosis of atherosclerosis (34). In gut inflammatory diseases, TNFα induced the expression of MMP1, which remodeled ECM to cause colon damage. MMP1 could also increase the expression of TNFα in a positive feedback mechanism to result in excessive mucosal damage (35). D’Haens et al. have developed 13 blood biomarkers to improve diagnosis of patients with Crohn’s disease, which included MMP1, MMP3, MMP9 (36). In joint inflammatory disease, synoviocytes expressed MMP1, MMP3, and MMP10 to increase their invasiveness to destruct the cartilage in rheumatoid arthritis (37). MMPs were involved in the proteolytic cleavage of important tissue components surrounding the joint and collagen fragments degradation, exposing various immunodominant epitopes for immune response (38). Besides, MMPs participate in cancer progression, metastasis and angiogenesis (9).
There is an increasing focus on the role of MMPs in renal diseases. Serum MMP1 was found to be significantly elevated in polycystic kidney disease (39). MMP3 could mediate shedding of kidney injury molecule-1 in renal tubular epithelial cells during acute kidney injury (40). The upregulation of MMP9 induced endothelial mesenchymal transition of tubular cells to promote kidney fibrosis (41). The increased expression of MMP10 was found specifically in the podocytes of injured glomeruli, further causing proteinuria and glomerulosclerosis (42). However, the relationship between MMPs and urolithiasis has been little investigated. Wu et al. found a high calcium concentration activated the ROS/NF-κB/MMP9 axis in tubular cells to facilitate osteoblastic transformation and crystals deposition (43). Also, MMP9 gene polymorphisms was also found to be related with urolithiasis (44).
We divided 29 RPs samples into two clusters based on the expression of hub DEMMPs. Cluster A exhibited high expression of hub DEMMPs while cluster B exhibited low expression. GO analysis showed DEGs were enriched in leukocyte mediated immunity, leukocyte chemotaxis, myeloid leukocyte migration, suggesting MMPs might participate in RPs progression through immunity and inflammation. KEGG analysis showed DEGs were enriched in inflammatory diseases, such as rheumatoid arthritis, inflammatory bowel disease, and these diseases were confirmed to be related with MMPs. The results of GSEA were also consistent with the above findings, as cluster A was enriched with activation of immune response, acute inflammatory response and adaptive immune response. Thus, we further conducted immune infiltration analysis and found that M1 macrophages were more abundant in cluster A. Markers of M1 macrophages, including NOS2, TLR2, and TLR4, were also higher in cluster A. All hub DEMMPs showed positive correlation tendency with M1 macrophages, although MMP9 and MMP10 did not reach statistical significance. Studies found that mice with a deficiency of MMP3 or MMP9 showed reduced macrophage infiltration in atherosclerosis in atherosclerotic plaques (45, 46). Similarly, our results indicated that increased expression of MMPs was accompanied by more M1 macrophages infiltration.
Macrophage function in CaOx stone formation has been a hot research topic in recent years (13). M0 macrophages can be polarized into two major phenotypes, pro-inflammatory M1 macrophages and anti-inflammatory M2 macrophages based on different conditions. The renal papillary tissues of stone formers exhibited higher expression of M1-related genes and lower expression of M2-related genes (47). Studies have confirmed that M1 macrophages could release pro-inflammatory molecules and induce renal damage and inflammation to promote crystals deposition and stone formation (47–49). Conversely, the M2 macrophages exert protective effects against stone formation owing to their potent capacity to phagocytose crystals (47, 50, 51). Hence, regulation of macrophage phenotypic transformation from M1 to M2 might be therapeutic targets.
Various studies have been conducted to elucidate the relationship between MMPs and macrophages. On the one hand, macrophages are an important source of MMPs, such as MMP9, MMP10, and MMP12. MMP9 from macrophages could promoted elastin degradation and accelerated disruption of atherosclerotic plaques (52). M1 macrophages could release MMP10 to induce pulmonary artery smooth muscle cells proliferation and migration, further leading to vascular remodeling and pulmonary arterial hypertension (53). MMP12, also called macrophage elastase, could mediate a variety of pathological processes. On the other hand, MMPs secreted by other cells can induce macrophage infiltration. For example, MMP9 induced macrophages recruitment through cleavage of SPP1 (54). Several MMPs, including MMP1, MMP3, MMP9, and MMP12, can process proTNF on macrophages to its active form, which might evoke the constitutive release of TNF from macrophages to induce tissue damage (55). We assumed that cluster A might recruit more M1 macrophages around RPs to exacerbate inflammatory responses and promote RPs progression and stone formation.
This study has several limitations. First, our results are mostly based on bioinformatics analysis and the underlying mechanisms behind the role of MMPs in urolithiasis remains unclear, which need to be further verified and clarified by in vitro and in vivo experiments. second, the sample size of the public dataset was still relatively small, which require further large samples validation.
Conclusion
In summary, we identified five hub DEMMPs implicated in urolithiasis via WGCNA and machine learning algorithms. What’s more, RPs samples were classified into two clusters based on their expression profile. Cluster A which exhibited higher expression of hub DEMMPs is enriched with immune response and inflammation through functional analysis. Besides, cluster A showed higher immune infiltration levels, especially increased infiltration of M1 macrophages and enhanced levels of inflammation compared to cluster B. Our findings offer a novel perspective on the involvement of MMPs in the development of RPs and stone formation for the first time, and provided potential targets for treatment and prevention.
Data availability statement
The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding authors.
Author contributions
S-YH, S-GW and B-LQ designed the study. S-YH and H-CJ performed the analysis. S-YH, H-CJ and W-CX wrote the manuscript. W-CX and H-SZ contributed to preparing the figures and tables. H-SZ and S-GW revised the manuscript. B-LQ provided the funding and supervised the study. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the National Natural Science Foundation of China (grant number 81900646).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Abbreviations
CaOx, Calcium oxalate; CaP, Calcium phosphate; RPs, Randall’s plaques; MMPs, Matrix metalloproteinases; ECM, Extracellular matrix; DEMMPs, Differentially expressed MMPs; WGCNA, Weighted gene co-expression network analysis; DEGs, Differentially expressed genes; GSEA, Gene set enrichment analysis; PPI, Protein-protein interaction; LASSO, Least absolute shrinkage and selection operator; SVM-RFE, Support vector machine recursive feature elimination; RF, Random forest; ROC, Receiver operating characteristic; AUC, Area under the ROC curve; CDF, Cumulative distribution function; PCA, Principal component analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ssGSEA, Single sample gene set enrichment analysis.
References
1. Sorokin I, Mamoulakis C, Miyazawa K, Rodgers A, Talati J, Lotan Y. Epidemiology of stone disease across the world. World J Urol (2017) 35(9):1301–20. doi: 10.1007/s00345-017-2008-6
2. D'Costa MR, Haley WE, Mara KC, Enders FT, Vrtiska TJ, Pais VM, et al. Symptomatic and radiographic manifestations of kidney stone recurrence and their prediction by risk factors: A prospective cohort study. J Am Soc Nephrol (2019) 30(7):1251–60. doi: 10.1681/asn.2018121241
3. Singh P, Enders FT, Vaughan LE, Bergstralh EJ, Knoedler JJ, Krambeck AE, et al. Stone composition among first-time symptomatic kidney stone formers in the community. Mayo Clin Proc (2015) 90(10):1356–65. doi: 10.1016/j.mayocp.2015.07.016
4. Khan SR, Pearle MS, Robertson WG, Gambaro G, Canales BK, Doizi S, et al. Kidney stones. Nat Rev Dis Primers (2016) 2:16008. doi: 10.1038/nrdp.2016.8
5. Daudon M, Bazin D, Letavernier E. Randall's plaque as the origin of calcium oxalate kidney stones. Urolithiasis (2015) 43 Suppl 1:5–11. doi: 10.1007/s00240-014-0703-y
6. Khan SR, Canales BK, Dominguez-Gutierrez PR. Randall's plaque and calcium oxalate stone formation: Role for immunity and inflammation. Nat Rev Nephrol (2021) 17(6):417–33. doi: 10.1038/s41581-020-00392-1
7. Khan SR, Canales BK. Unified theory on the pathogenesis of randall's plaques and plugs. Urolithiasis (2015) 43 Suppl 1(0 1):109–23. doi: 10.1007/s00240-014-0705-9
8. Nagase H, Visse R, Murphy G. Structure and function of matrix metalloproteinases and TIMPs. Cardiovasc Res (2006) 69(3):562–73. doi: 10.1016/j.cardiores.2005.12.002
9. de Almeida LGN, Thode H, Eslambolchi Y, Chopra S, Young D, Gill S, et al. Matrix metalloproteinases: From molecular mechanisms to physiology, pathophysiology, and pharmacology. Pharmacol Rev (2022) 74(3):712–68. doi: 10.1124/pharmrev.121.000349
10. Laronha H, Caldeira J. Structure and function of human matrix metalloproteinases. Cells (2020) 9(5):1076. doi: 10.3390/cells9051076
11. Zakiyanov O, Kalousová M, Zima T, Tesař V. Matrix metalloproteinases and tissue inhibitors of matrix metalloproteinases in kidney disease. Adv Clin Chem (2021) 105:141–212. doi: 10.1016/bs.acc.2021.02.003
12. Taguchi K, Hamamoto S, Okada A, Unno R, Kamisawa H, Naiki T, et al. Genome-wide gene expression profiling of randall's plaques in calcium oxalate stone formers. J Am Soc Nephrol (2017) 28(1):333–47. doi: 10.1681/asn.2015111271
13. Taguchi K, Okada A, Unno R, Hamamoto S, Yasui T. Macrophage function in calcium oxalate kidney stone formation: A systematic review of literature. Front Immunol (2021) 12:673690. doi: 10.3389/fimmu.2021.673690
14. Parks WC, Wilson CL, López-Boado YS. Matrix metalloproteinases as modulators of inflammation and innate immunity. Nat Rev Immunol (2004) 4(8):617–29. doi: 10.1038/nri1418
15. Langfelder P, Horvath S. WGCNA: An r package for weighted correlation network analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559
16. Yip AM, Horvath S. Gene network interconnectedness and the generalized topological overlap measure. BMC Bioinf (2007) 8:22. doi: 10.1186/1471-2105-8-22
17. Ravasz E, Somera AL, Mongru DA, Oltvai ZN, Barabási AL. Hierarchical organization of modularity in metabolic networks. Science (2002) 297(5586):1551–5. doi: 10.1126/science.1073374
18. Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Society. Ser B (Methodological) (1996) 58(1):267–88. doi: 10.1111/j.2517-6161.1996.tb02080.x
19. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw (2010) 33(1):1–22. doi: 10.18637/jss.v033.i01
20. Sanz H, Valim C, Vegas E, Oller JM, Reverter F. SVM-RFE: Selection and visualization of the most relevant features through non-linear kernels. BMC Bioinf (2018) 19(1):432. doi: 10.1186/s12859-018-2451-4
21. Karatzoglou A, Smola A, Hornik K, Zeileis A. Kernlab - an S4 package for kernel methods in r. J Stat Software (2004) 11(9):1–20. doi: 10.18637/jss.v011.i09
24. Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, et al. pROC: an open-source package for r and s+ to analyze and compare ROC curves. BMC Bioinf (2011) 12:77. doi: 10.1186/1471-2105-12-77
25. Wilkerson MD, Hayes DN. ConsensusClusterPlus: A class discovery tool with confidence assessments and item tracking. Bioinformatics (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170
26. Yu G, Wang LG, Han Y, He QY. clusterProfiler: An r package for comparing biological themes among gene clusters. Omics (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118
27. Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell (2015) 160(1-2):48–61. doi: 10.1016/j.cell.2014.12.033
28. Hänzelmann S, Castelo R, Guinney J. GSVA: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7
29. Randall A. The origin and growth of renal calculi. Ann Surg (1937) 105(6):1009–27. doi: 10.1097/00000658-193706000-00014
30. Kim HE, Dalal SS, Young E, Legato MJ, Weisfeldt ML, D'Armiento J. Disruption of the myocardial extracellular matrix leads to cardiac dysfunction. J Clin Invest (2000) 106(7):857–66. doi: 10.1172/jci8040
31. Hadi T, Boytard L, Silvestro M, Alebrahim D, Jacob S, Feinstein J, et al. Macrophage-derived netrin-1 promotes abdominal aortic aneurysm formation by activating MMP3 in vascular smooth muscle cells. Nat Commun (2018) 9(1):5022. doi: 10.1038/s41467-018-07495-1
32. Dwivedi A, Slater SC, George SJ. MMP-9 and -12 cause n-cadherin shedding and thereby beta-catenin signalling and vascular smooth muscle cell proliferation. Cardiovasc Res (2009) 81(1):178–86. doi: 10.1093/cvr/cvn278
33. Matilla L, Roncal C, Ibarrola J, Arrieta V, García-Peña A, Fernández-Celis A, et al. A role for MMP-10 (Matrix metalloproteinase-10) in calcific aortic valve stenosis. Arterioscler Thromb Vasc Biol (2020) 40(5):1370–82. doi: 10.1161/atvbaha.120.314143
34. Hu W, Wei R, Wang L, Lu J, Liu H, Zhang W. Correlations of MMP-1, MMP-3, and MMP-12 with the degree of atherosclerosis, plaque stability and cardiovascular and cerebrovascular events. Exp Ther Med (2018) 15(2):1994–98. doi: 10.3892/etm.2017.5623
35. Wang YD, Mao JW. Expression of matrix metalloproteinase-1 and tumor necrosis factor-alpha in ulcerative colitis. World J Gastroenterol (2007) 13(44):5926–32. doi: 10.3748/wjg.v13.i44.5926
36. D'Haens G, Kelly O, Battat R, Silverberg MS, Laharie D, Louis E, et al. Development and validation of a test to monitor endoscopic activity in patients with crohn's disease based on serum levels of proteins. Gastroenterology (2020) 158(3):515–26.e10. doi: 10.1053/j.gastro.2019.10.034
37. Tolboom TC, Pieterman E, van der Laan WH, Toes RE, Huidekoper AL, Nelissen RG, et al. Invasive properties of fibroblast-like synoviocytes: Correlation with growth characteristics and expression of MMP-1, MMP-3, and MMP-10. Ann Rheum Dis (2002) 61(11):975–80. doi: 10.1136/ard.61.11.975
38. Dufour A, Overall CM. Missing the target: Matrix metalloproteinase antitargets in inflammation and cancer. Trends Pharmacol Sci (2013) 34(4):233–42. doi: 10.1016/j.tips.2013.02.004
39. Nakamura T, Ushiyama C, Suzuki S, Ebihara I, Shimada N, Koide H. Elevation of serum levels of metalloproteinase-1, tissue inhibitor of metalloproteinase-1 and type IV collagen, and plasma levels of metalloproteinase-9 in polycystic kidney disease. Am J Nephrol (2000) 20(1):32–6. doi: 10.1159/000013552
40. Lim AI, Chan LY, Lai KN, Tang SC, Chow CW, Lam MF, et al. Distinct role of matrix metalloproteinase-3 in kidney injury molecule-1 shedding by kidney proximal tubular epithelial cells. Int J Biochem Cell Biol (2012) 44(6):1040–50. doi: 10.1016/j.biocel.2012.03.015
41. Strutz F, Zeisberg M, Ziyadeh FN, Yang CQ, Kalluri R, Müller GA, et al. Role of basic fibroblast growth factor-2 in epithelial-mesenchymal transformation. Kidney Int (2002) 61(5):1714–28. doi: 10.1046/j.1523-1755.2002.00333.x
42. Zuo Y, Wang C, Sun X, Hu C, Liu J, Hong X, et al. Identification of matrix metalloproteinase-10 as a key mediator of podocyte injury and proteinuria. Kidney Int (2021) 100(4):837–49. doi: 10.1016/j.kint.2021.05.035
43. Wu Y, Zhang J, Li C, Hu H, Qin B, Wang T, et al. The activation of ROS/NF-κB/MMP-9 pathway promotes calcium-induced kidney crystal deposition. Oxid Med Cell Longev (2021) 2021:8836355. doi: 10.1155/2021/8836355
44. Mehde AA, Mehdi WA, Yusof F, Raus RA, Abidin ZAZ, Ghazali H, et al. Association of MMP-9 gene polymorphisms with nephrolithiasis patients. J Clin Lab Anal (2018) 32(1):e22173. doi: 10.1002/jcla.22173
45. Silence J, Lupu F, Collen D, Lijnen HR. Persistence of atherosclerotic plaque but reduced aneurysm formation in mice with stromelysin-1 (MMP-3) gene inactivation. Arterioscler Thromb Vasc Biol (2001) 21(9):1440–5. doi: 10.1161/hq0901.097004
46. Luttun A, Lutgens E, Manderveld A, Maris K, Collen D, Carmeliet P, et al. Loss of matrix metalloproteinase-9 or matrix metalloproteinase-12 protects apolipoprotein e-deficient mice against atherosclerotic media destruction but differentially affects plaque growth. Circulation (2004) 109(11):1408–14. doi: 10.1161/01.Cir.0000121728.14930.De
47. Taguchi K, Okada A, Hamamoto S, Unno R, Moritoki Y, Ando R, et al. M1/M2-macrophage phenotypes regulate renal calcium oxalate crystal development. Sci Rep (2016) 6:35167. doi: 10.1038/srep35167
48. Taguchi K, Okada A, Hamamoto S, Iwatsuki S, Naiki T, Ando R, et al. Proinflammatory and metabolic changes facilitate renal crystal deposition in an obese mouse model of metabolic syndrome. J Urol (2015) 194(6):1787–96. doi: 10.1016/j.juro.2015.07.083
49. Dominguez-Gutierrez PR, Kusmartsev S, Canales BK, Khan SR. Calcium oxalate differentiates human monocytes into inflammatory M1 macrophages. Front Immunol (2018) 9:1863. doi: 10.3389/fimmu.2018.01863
50. Taguchi K, Okada A, Kitamura H, Yasui T, Naiki T, Hamamoto S, et al. Colony-stimulating factor-1 signaling suppresses renal crystal formation. J Am Soc Nephrol (2014) 25(8):1680–97. doi: 10.1681/asn.2013060675
51. Anders HJ, Suarez-Alvarez B, Grigorescu M, Foresto-Neto O, Steiger S, Desai J, et al. The macrophage phenotype and inflammasome component NLRP3 contributes to nephrocalcinosis-related chronic kidney disease independent from IL-1-mediated tissue injury. Kidney Int (2018) 93(3):656–69. doi: 10.1016/j.kint.2017.09.022
52. Gough PJ, Gomez IG, Wille PT, Raines EW. Macrophage expression of active MMP-9 induces acute plaque disruption in apoE-deficient mice. J Clin Invest (2006) 116(1):59–69. doi: 10.1172/jci25074
53. Chi PL, Cheng CC, Hung CC, Wang MT, Liu HY, Ke MW, et al. MMP-10 from M1 macrophages promotes pulmonary vascular remodeling and pulmonary arterial hypertension. Int J Biol Sci (2022) 18(1):331–48. doi: 10.7150/ijbs.66472
54. Tan TK, Zheng G, Hsu TT, Lee SR, Zhang J, Zhao Y, et al. Matrix metalloproteinase-9 of tubular and macrophage origin contributes to the pathogenesis of renal fibrosis via macrophage recruitment through osteopontin cleavage. Lab Invest (2013) 93(4):434–49. doi: 10.1038/labinvest.2013.3
Keywords: urolithiasis, bioinformatics, MMPs, immune infiltration, macrophages
Citation: Hong S-Y, Jiang H-C, Xu W-C, Zeng H-S, Wang S-G and Qin B-L (2023) Bioinformatics analysis reveals the potential role of matrix metalloproteinases in immunity and urolithiasis. Front. Immunol. 14:1158379. doi: 10.3389/fimmu.2023.1158379
Received: 03 February 2023; Accepted: 06 March 2023;
Published: 15 March 2023.
Edited by:
Guohua Zeng, First Affiliated Hospital of Guangzhou Medical University, ChinaReviewed by:
Zhong Wu, Fudan University, ChinaKehua Jiang, Guizhou Provincial People’s Hospital, China
Copyright © 2023 Hong, Jiang, Xu, Zeng, Wang and Qin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Shao-Gang Wang, c2d3YW5ndGptQDE2My5jb20=; Bao-Long Qin, cWluYmFvbG9uZ0BodXN0LmVkdS5jbg==