The mRNA microarray dataset GSE73680 (Taguchi et al., 2017) and lncRNA microarray dataset GSE117518 (Zhu et al., 2021) were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE73680 dataset included 24 RPs from patients with CaOx stones and six controls with normal papillary tissues. The GSE117518 dataset included three RPs from patients with CaOx stones and three controls with normal papillary tissues. The details of both datasets are presented in Table 1. Probe names were transformed into gene symbols according to platform annotation information. Moreover, immune-related genes were obtained from the Immunology Database and Analysis Portal (IMMPORT) database (http://www.immport.org/) (Bhattacharya et al., 2014).

Subsequently, DE-mRNAs and lncRNAs were analyzed and compared between RPs and normal–papillary tissue controls using the “limma” package (Ritchie et al., 2015) in the R software (http://www.r-project.org). mRNAs that met the criteria of |log₂FC|>1 and p < 0.01 were considered as DE-mRNAs, and lncRNAs that met the criteria of |log₂FC|>0.58 and p < 0.01 were considered as DE-lncRNA. The “ggplot2” package was used to draw heatmaps and volcano plots for data visualization.

Potential DE-lncRNA-targeted miRNAs were predicted using the miRcode database (http://mircode.org/) (Jeggari et al., 2012). Subsequently, miRNA–mRNA pairs were analyzed using TargetScan (http://www.targetscan.org/vert_72/) (Agarwal et al., 2015), miRTarBase (https://mirtarbase.cuhk.edu.cn/php/index.php) (Huang et al., 2020), and miRDB (http://mirdb.org/) (Chen and Wang, 2020) databases. mRNAs that were found in at least two databases were considered as candidate targets of miRNAs.

The Venn method was used to analyze overlapping genes. Intersections between DE-mRNA and DE-lncRNA-targeted mRNAs, as well as lncRNA-related DE-mRNAs and immune-related genes were identified using the Venny version 2.1 online tool (https://bioinfogp.cnb.csic.es/tools/venny/index.html).

To explore the functions of lncRNA-related DE-mRNAs, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted using the “org.Hs.eg.db” and “ClusterProfiler” packages (Yu et al., 2012) in the R software. An adjusted p < 0.05 was considered statistically significant. Subsequently, the STRING database (https://string-db.org/) (Szklarczyk et al., 2019) was used to determine the relationship between the DE-mRNAs, and Cytoscape software (https://cytoscape.org) was used to develop the PPI network.

After identifying immune-related and lncRNA-related DE-mRNAs, the interaction between lncRNAs, miRNAs, and mRNAs was confirmed as described in item 3.2. Subsequently, the immune-related lncRNA–miRNA–mRNA ceRNA network was developed using the R software. The “ggalluvial” package was used to draw a sankey diagram for data visualization.

To estimate the abundance of 22 types of immune cell types in Randall’s plaques and normal–papillary tissue controls, the mRNA microarray dataset GSE73680 was uploaded to the platform of CIBERSORT (http://cibersort.stanford.edu/) (Newman et al., 2015). Only samples that had a CIBERSORT algorithm output of p < 0.05 were considered for further analysis. Histograms and heatmaps were drawn to show the rate of immune cell infiltration in different samples. Co-expression patterns in immune-related DE-mRNAs and infiltrating immune cells were analyzed using Pearson’s correlation coefficient. Subsequently, the Wilcoxon rank-sum test was performed to compare differentially infiltrating immune cells between RPs in patients with CaOx stones and controls with normal papillary tissues. The relationship between DE-mRNA expression and the fractions of macrophages and mast cells was also investigated using the Wilcoxon test. Results were visualized using the “heatmap” and “vioplot” packages in the R software.

Thirty male C57BL/6J mice weighing 22–25 g and aging 6–8 weeks were acquired from the Center of Experimental Animals at the Renmin hospital of Wuhan University, Hubei, China. The mice were acclimatized in the animal house of our institution at a steady temperature of 22 ± 2°C and humidity of 40–70% on a 12/12-h light–dark cycle and with free access to water and feed. The animal experiments were conducted according to the Guide for the Care and Use of Laboratory Animals, and the study protocol was approved by the Laboratory Animal Welfare and Ethics Committee of the Renmin hospital of Wuhan University (approval number: WDRM-20200604).

According to previous publications (Okada et al., 2007; Usami et al., 2018), the mice were intraperitoneally injected with 80 or 120 mg/kg of glyoxylate (Sigma–Aldrich; St. Louis, MO, United States) daily for 6 or 12 days to establish a CaOx kidney stone model. Mice were randomly assigned to the five following dosage groups (n = 6): control, 80 mg/kg of glyoxylate for 6 days, 120 of mg/kg glyoxylate for 6 days, 80 mg/kg of glyoxylate for 12 days, and 120 mg/kg of glyoxylate for 12 days groups. After 6 or 12 days, the mice were sacrificed, and kidneys were removed for analyses.

Human renal tubular epithelial cell line (HK-2) cells were provided by Stem Cell Bank, Chinese Academy of Sciences, Shanghai, China. HK-2 cells were cultured in an MEM medium supplemented with 10% fetal bovine serum (Gibco, Waltham, MA, United States) and 1% antibiotics (penicillin/streptomycin). The cells were maintained at 37°C under a humidified atmosphere with 5% CO₂. Oxalate was purchased from Sigma–Aldrich and dissolved in the culture medium. Subsequently, the cells were cultured in six-well plates, and 0.25 mM or 0.5 mM of oxalate were added for 48 h.

After fixation in 4% paraformaldehyde, kidneys were imbedded in paraffin and were cut into 5-μm slices. HE staining was performed to assess the histopathological kidney tubular injuries as previously described (Dong et al., 2019). Injuries were scored as follows: 0, no tubular injury; 1, <10% tubular damage; 2, 10–25% tubular damage; 3, 25–50% tubular damage; 4, 50–74% tubular damage; and 5, >75% tubular damage. Subsequently, crystals were detected using Von Kossal staining, as previously described (Wang et al., 2019). The crystal deposition area was quantified using Image J software. Renal fibrosis was verified using Masson trichrome staining, and the collagen fiber deposition area on kidney sections was quantified using Image J software.

The protein expression levels of TGF-β, iNOS, and CD206 were analyzed using immunohistochemical and immunofluorescence staining. Antibodies (i.e., TGF-β [21898-1-AP], iNOS [18985-1-AP], and CD206 [60143-1-Ig]) were purchased from Proteintech (Chicago, IL, United States). All procedures were conducted according to the recommendations of the manufacturer. By comparing the positive area between groups using microscopy, figures were analyzed using Image J software.

Mast cells were detected using Toluidine blue staining as previously described (Zhang et al., 2017). Mast cells were identified using purple granules, and activated mast cells were characterized by disgorged and loosely packed granules. Activated mast cells per field were counted at a magnification of 400 ×.

Total RNA was extracted from HK-2 cells using TRIzol reagent (Invitrogen Life Technologies, Carlsbad, CA, United States), and RNA purity was measured using spectrophotometry. RNAs were reverse transcribed into cDNAs using the Takara RNA PCR kit (Takara Biotechnology, Shiga, Japan) according to the instructions of the manufacturer. Subsequently, cDNA was amplified by RT-qPCR using an Applied Biosystems SYBR Green mix kit (Applied Biosystems, Foster City, CA, United States). GAPDH was used as an internal reference for lncRNAs and mRNAs, while U6 was used as a reference for miRNA. The primers used for these reactions are shown in Supplementary Table S1. The reactions were measured on the ABI 7900 Real-Time PCR system (Applied Biosystems Life Technologies), and the 2^−ΔΔCT method was used for analysis.

All data are presented as the mean ± SD. Statistical analysis was conducted using SPSS version 19.0 (SPSS Inc., Chicago, IL, United States). Student’s t-test was used to compare differences between groups. A p-value of <0.05 was considered statistically significant. All experiments were performed at least three times.

To clarify the process of this research, a schematic representation is presented in Figure 1. Original data were downloaded from the GSE73680 and GSE117518 datasets in the GEO database. In the GSE73680 dataset, RNA-seq data of 24 RPs from patients with CaOx stones and from six controls with normal papillary tissues were analyzed using criteria of |log₂FC|>1 and p < 0.01. A total of 2,340 DE-mRNAs (2,098 upregulated and 242 downregulated) were compared between RPs and normal papillary tissues. In the GSE117518 dataset, the RNA-seq data of three RPs from patients with CaOx stones and three normal papillary tissues were analyzed using criteria of |log₂FC|>0.58 and p < 0.01. A total of 929 DE-lncRNAs (587 upregulated and 342 downregulated) were identified. Corresponding heatmaps and volcano plots are shown in Figure 2. Details of datasets are presented in Table 1.

To establish the ceRNA network, DE-lncRNAs were further analyzed. Potential DE-lncRNA-targeted miRNAs were predicted using the miRcode database. Subsequently, miRNA–mRNA pairs were analyzed using the TargetScan, miRTarBase, and miRDB databases. A total of 197 miRNAs and 8,457 mRNAs were predicted. Subsequently, the Venn method was used to analyze the intersection between DE-mRNA and DE-lncRNA-targeted mRNAs (Figure 3A). Consequently, 278 overlapping lncRNA-related DE-mRNAs were identified. To determine the functions of lncRNA-related DE-mRNAs, GO and KEGG enrichment analyses were conducted (Figures 3B,C). A biological process analysis showed that lncRNA-related DE-mRNAs were significantly enriched in extracellular matrix organization, cellular calcium ion homeostasis, and regulation of cellular response to growth factor stimulus. A cellular component analysis showed that lncRNA-related DE-mRNAs were mostly enriched in collagen-containing extracellular matrix and endoplasmic reticulum lumen. A molecular function (MF) analysis showed that lncRNA-related DE-mRNAs were mostly enriched in channel activity and extracellular matrix structural constituent. The KEGG pathway enrichment analysis showed that lncRNA-related DE-mRNAs were significantly enriched in PI3K-Akt signaling pathway, focal adhesion, and extracellular matrix-receptor interaction. Details of GO and KEGG enrichment analyses are presented in Tables 2, 3. The PPI network of lncRNA-related DE-mRNAs is shown in Supplementary Figure S1.

To construct the immune-related ceRNA network, the Venn method was used to analyze the intersection between lncRNA-related DE-mRNAs and immune-related genes obtained from the IMMPORT database. Consequently, 20 overlapping immune-related DE-mRNAs (12 upregulated and eight downregulated) were identified (Figure 4A). A co-expression analysis of immune-related DE-mRNAs was performed (Figure 4B). Subsequently, immune-related DE-mRNAs and their paired miRNAs and lncRNAs were chosen to develop the ceRNA regulatory network (Figure 4C). In total, the immune-related ceRNA network contained 10 lncRNAs, 23 miRNAs, and 20 mRNAs.

The composition of 22 infiltrating immune cells in RPs in patients with CaOx stones and controls with normal papillary tissues were estimated using the CIBERSORT algorithm (Figures 5A,B). The relationships among these 22 immune cells are presented in Figure 5C. M1 macrophages were positively correlated with resting dendritic cells (R = 0.70). M2 macrophages were positively correlated with eosinophils (R = 0.52). Activated mast cells activated were positively correlated with neutrophils (R = 0.59). Resting mast cells were positively correlated with activated NK cells (R = 0.55) and negatively correlated with resting dendritic cells (R = –0.45). The differential proportion of infiltrating immune cells between RPs in patients with CaOx stones and in controls with normal papillary tissues was analyzed. As shown in Figure 5D, compared with the RPs in controls, M2 macrophages (p = 0.038) and resting mast cells (p = 0.019) were significantly downregulated and M1 macrophages (p = 0.49) and activated mast cells (p = 0.296) were significantly upregulated in the RPs in patients with kidney stones.

For further analysis, DE-mRNAs were divided to the high expression and low expression groups. The correlation between infiltrating immune cells and DE-mRNAs expression was estimated using the Wilcoxon test, and significantly correlated pairs with p < 0.05 are shown in Figure 6. Results indicated that the expression of IL-13, OGN, and VEGFC was significantly negatively correlated with the proportion of M1 macrophages (p = 0.011, p = 0.002, and p = 0.05, respectively), whereas the expression of VAV2 was significantly positively correlated with proportion of M1 macrophages (p = 0.038). The expression of ADM2, CCL7, FGF18, FGF21, CCR9, LEP, ROBO2, and VAV2 was significantly negatively correlated with the proportion of resting mast cells (p = 0.011, p = 0.049, p = 0.001, p = 0.002, p < 0.001, p = 0.016, p = 0.007, p = 0.02, and p = 0.002, respectively).

To validate the aforementioned pathway and differentially infiltrating immune cells, kidney tubular injury, crystal deposition, fibrosis level, and macrophage and mast cell infiltration were evaluated in a glyoxylate (Gly)-induced hyperoxaluric mouse model. As shown in Figure 7, kidney tubular injury and crystal deposition were aggravated as treatment concentration and time of glyoxylate increased. Tubular injury and crystals were markedly worse in the day–12 Gly-treated groups than in the day–6 Gly-treated groups. Moreover, tubular injury and crystals were markedly worse in the 120–mg/kg Gly-treated mice than in the 80–mg/kg Gly-treated mice in both day-6 and day-12 groups. Fibrosis and collagen fiber deposition were evaluated using Masson staining and immunohistochemical staining of TGF-β. Results have shown that collagen fiber depositions and TGF-β-positive areas were significantly more frequent in the Gly-treated groups than in the control group in a dose- and time-dependent manner; these results are consistent with those shown in Figure 3B.

Subsequently, the immunofluorescence staining of macrophage-related molecules iNOS (M1) and CD206 (M2) showed that M1 macrophage infiltration significantly increased as kidney stones aggravated, whereas M2 macrophage infiltration increased in the early stages (day 6) and decreased as kidney stones progressed (day 12). Toluidine blue staining showed that activated mast cell infiltration significantly increased in the kidneys of mice with stone formation. As treatment concentration and time of Gly increased, activated mast cells concomitantly increased. These immune cell infiltration results are consistent with our findings shown in Figure 5.

Through literature review and co-expression analysis of infiltrating immune cells, we established the immune-related hub ceRNA network, comprising 2 lncRNAs, 3 miRNAs, and 2 mRNAs (Figure 8A). Details of the immune-related hub ceRNA network developed from the GSE73680 and GSE117518 datasets are presented in Table 4. To validate the immune-related hub ceRNA network in kidney stones, RT-qPCR was used to detect the expression levels of the hub genes. As shown in Figure 8B, treatment with 0.25 and 0.5 mM oxalate for 48 h significantly upregulated the expression levels of NEAT1, PVT1, CCL7, and ROBO2 but downregulated the expression levels of hsa-miR-23b-3p, hsa-miR-429, and hsa-miR-139-5p in the HK-2 cell line in a dose-dependent manner. These results are consistent with the findings of GEO datasets.