MME and PTPRC: key renal biomarkers in lupus nephritis

Data sources

For the current study, the datasets, namely GSE200306 (Parikh et al., 2022) (Platforms: GPL21847 (HG-U133_Plus_2)) have been retrieved from Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/). This dataset has been meticulously curated, encompassing both clinical and histological data for all patients. Initially, the GSE200306 dataset included 58 paired kidney biopsies including both first biopsy and second biopsy from individuals diagnosed with proliferative lupus nephritis. These biopsies were classified as either class III, class IV, or class III/IV+V, in which glomeruli and tubulointerstitium was isolated separately. Our study specifically analyzed only the first biopsy data of LN from this dataset including 34 LN glomerular samples and 45 LN tubulointerstitial samples. Additionally, nine glomerular samples and 10 tubulointerstitium samples from pre-implantation donor kidney biopsies were used as healthy controls (HC). Overall, three clusters from the datasets were delineated for bioinformatic analysis: glomerular cluster (LN_G group, 34 LN and 9 HC), tubulointerstitium cluster pairs (LN_T group, 45 LN and 10 HC), and integrated total renal cluster (LN group, 79 LN and 19 HC).

Identification of DEGs

The GEO dataset (GSE200306) was downloaded using the “GEOquery” package in R software. Expression data and phenotypic data were extracted from the dataset. First, we utilized the normalizeBetweenArrays function from the “limma” package (Ritchie et al., 2015) to normalize the microarray data ensuring that the expression values are comparable across different arrays. Then, a design matrix was created to represent the different groups under comparison. This matrix was essential for fitting the linear model. we analyzed differentially expressed genes (DEGs) in three distinct clusters: total kidney (LN group), glomeruli (LN_G group), and tubulointerstitium (LN_T group) from primary LN patients and healthy controls. DEGs were identified based on the following criteria: adjust. P < 0.05 and fold change (FC) >1. The “ggplot2” package was employed to visualize the differential DEGs results including volcano plot, heatmap and boxplot. The full code for R can be found in the R code supplemental file.

Functional and pathway enrichment analysis

We employed the “clusterProfiler” package in R for conducting functional enrichment analysis of the differentially expressed genes (DEGs) based on Gene Ontology (GO, http://www.geneontology.org) and the Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.kegg.jp). In the context of the GO analysis, we identified and explored three distinct categories: biological process (BP), cellular component (CC), and molecular function (MF) (The Gene Ontology C, 2019), shedding light on the biological processes associated with these DEGs. We also explored potential signaling pathways using the KEGG analysis.

Vene plots and PPI networks construction

To investigate the core differentially expressed genes (DEGs) in primary LN across different renal compartments, we analyzed the total renal group (LN group), the glomeruli parts (LN_G group), and the renal tubulointerstitium (LN_T group). Using an online tool (https://www.xiantaozi.com/), we generated Venn plots to visualize the overlaps and distinctions among these groups. To evaluate the relationships between the DEGs in the dataset, we established a protein-protein interaction (PPI) network using the Search Tool for the Retrieval of Interacting Genes (STRING v11.5) (Szklarczyk et al., 2021). Subsequently, we integrated the results obtained from the STRING database into Cytoscape software (v3.9.1) to visualize the protein-protein interactions (PPIs) within the statistically significant DEGs.

NephroSeq datasets verification

NephroSeq database (https://nephroseq.org/) is an free online platform designed for the exploration and analysis of gene expression data specifically related to nephrology. NephroSeq contains a vast collection of gene expression datasets from human nephropathy studies, such as diabetic nephropathy, IgA nephropathy, lupus nephritis. NephroSeq was used in our study to verify the changing trends and differences in the expression of seven DEGs identified from the GSE200306 dataset. By leveraging the extensive datasets available in NephroSeq, we were able to cross-validate our findings and gain additional insights into the behavior of these genes in LN datasets. A heat map displayed the p-values for the up and down expression of the seven DEGs, providing a clear visual representation of the statistical significance of gene expression changes.

ROC curve for screening the key biomarker

Using the outcome indicators in original sample data of GSE200306 to analyze the most important DEGs of glomeruli parts (LN_G group), and renal tubulointerstitium parts (LN_T group) respectively in primary LN by prism (version 8.0.0). And we also plotted ROC curves of common DEGs in the LN_G group and LN_T group to identify more diagnostically valuable common DEGs. When the value of Area Under Curve (AUC)>0.8, it means that the performance of the classifier is better, and the prediction result of the classifier is more reliable.

Animals experiment in MRL/lpr mice

Purchased 8-week-old SPF-level female MRL/lpr lupus mice (n = 3) and 8-week-old SPF-level female healthy control mice C57BL/6 (n = 3) from the SPF Biotechnology Co., Ltd. Mice were fed standard laboratory chow and supplied drinking water ad libitum in Department of Zoology, Central South University. All animal experiments were performed in accordance with the Guidelines for the Care and Use of Laboratory Animals published by the US National Institutes of Health (NIH). No animal intervention experimental study was carried out in this study. The criteria for euthanizing animals in our studies include: (1) Severe illness or injury that cannot be effectively treated; (2) significant and unrelievable pain or distress; (3) failure to thrive or meet basic welfare requirements; (4) development of unexpected complications that compromise the animal’s well-being; (5) completion of the experimental procedures or endpoints. If any of the aforementioned situations arise, we will euthanize the animals using carbon dioxide (CO2) euthanasia. Fill CO2 at an even rate to achieve a CO2 replacement rate of 30% to 70% of the container volume per minute, then animals lose consciousness and die within 2 to 3 min. Executed the mice when 25 weeks old and collected plasma, urine, and kidney samples for further analysis. Plasma samples were used to measure ds-DNA concentration, urine samples were analyzed for protein quantification, and kidney samples were used to assess mRNA and protein expression levels. Sampling procedures were conducted under light anesthesia using isoflurane to minimize discomfort and ensure the welfare of the animals. No animals survived at the conclusion of the experiment. The animal experiments mentioned above have been approved by Central South University’s Animal Welfare Ethics Review (Number: CSU-2022-0133).

Validation of key DEGs in MRL/lpr mice

The Bradford method was used to detect routine urine protein quantification. Plasma ds-DNA concentration was detected by mouse anti-double stranded DNA antibody (IgG) ELISA Kit (CSB-E11194m; CUSABIO). RT-qPCR was used to detect the mRNA expression in kidney tissue of five key DEGs screened in GSE200306. Primers of RT-qPCR were designed using primer 3 (https://primer3.ut.ee/), and the primer sequences are shown in Table S1. The reverse transcription was performed using HiScript II Q RT SuperMix for qPCR (Vazyme, R223), and the RT-qPCR process was carried out with ChamQ Universal SYBR qPCR Master Mix (Q711; Vazyme). The protein expression difference of DEGs in kidney was observed by immunohistochemistry, in which 10 different parts of each stained section were selected under a 40x optical microscope, and the average optical density was calculated using ImageJ2 software (version 2.3.0). Experimental antibodies were purchased from Servicebio. Statistical analysis was conducted using GraphPad Prism software (version 8.0.0). Unpaired t-test analysis was utilized for between-group comparisons, with statistical significance set at P < 0.05.

Identification of DEGs From the dataset

Our study utilized the gene expression profiles from the GSE200306 dataset obtained from the GEO database, which were originally submitted by Parikh et al. (2022), with platforms GPL21847 (HG-U133_Plus_2). After obtaining the original data, R software was used to control the quality of the data, which showed that the data of the LN group and the control group had good uniformity (Fig. 1A). A total of 497 genes were detected between the two groups, and the volcano plots showed significant differences in the minority (Fig. 1B). In the total renal of primary LN groups (LN group), we identified 17 hot DEGs (Fig. 1C). In comparison, 25 DEGs were discovered in the glomeruli group (LN_G group), and 24 DEGs were identified in the renal tubulointerstitium group (LN_T group) when comparing with control groups (Figs. 1D and 1E). These findings underscore the distinct molecular signatures present in different renal compartments.

Key DEGs and function analysis

Using Venn plot analysis, a total of seven significantly differential genes were shared among the three comparison groups (Fig. 2A). In these three groups of comparisons, the changes in the above seven DEGs are consistent, in which PTPRC and MX1 were the up-regulated genes, while EGR1, MME, RORC, ZBTB16, FKBP5 were the down-regulated DEGs. Table 1 showed the detail changes of the seven common DEGs in the total renal group (LN group). We used the STRING tool to discover the protein interactions among the common differential DEGs, as a result, PTPRC, MME and MX1 occupied the most central position in the network (Fig. 2B).

We performed functional enrichment analysis of DEGs between the LN and HC groups concentrating on KEGG and GO pathways. KEGG showed hematopoietic cell lineage, TH1/2/17 cell differentiation was the most related pathway in the LN pathogenesis (Fig. 2C). GO analysis demonstrated that T cell differentiation, lymphocyte differentiation and mononuclear cell differentiation participated in the key parts of the biological process of LN (Fig. 2D).

NephroSeq datasets analysis

Three datasets in the NephroSeq database (https://nephroseq.org/) from two published articles (Berthier et al., 2012; Peterson et al., 2004) and not yet published by ERCB in 2018/04/01 (see Data Availability section) were included in the analysis. The three databases contain both glomerular and tubulointerstitium gene clusters of lupus nephritis, in which total 62 subgroups from six Lupus datasets were analysed in our study, including Berthier Lupus Glom (n = 46), Berthier Lupus Mouse Kidney (n = 68), Berthier Lupus TubIntn (n = 47), ERCB Lupus TubInt (n = 41), ERCB Lupus Glom (n = 32) and Peterson Lupus Glom (n = 31) (Tables S2 and S3). In the comparison of 62 subgroups, the up-regulated gene (PTPRC, MX1) and the down-regulated DEGs (EGR1, MME, RORC, ZBTB16, FKBP5) were verified to have mostly consistent change trends in all Lupus Nephritis vs. Normal Kidneys group with the results of GSE200306 (Fig. 3). In addition, proteinuria status and LN pathological class also showed significantly different gene expression, especially in PTPRC and MME.

ROC curve for screening the key biomarker

The diagnostic ROC curves of DEGs were analyzed based on glomerulus and tubules respectively (Table S4). Among these, the seven common DEGs in the glomeruli group (LN_G) and the renal tubulointerstitium group (LN_T), which were identified in the Venn diagrams and verified in the NephroSeq database, were used to plot ROC curves and identify the diagnostic renal biomarkers for LN (Fig. 4). Consequently, five common DEGs including PTPRC, MX1, EGR1, MME and RORC, were found to be more promising biomarkers in both glomerulus (Fig. 4A) and tubulointerstitium (Fig. 4B) with AUC values greater than 0.80.

Key renal biomarkers verified by MRL/lpr mice

Using SPF-level female MRL/lpr lupus mice (n = 3) and female C57BL/6 (n = 3) to make a verification for the five key DEGs screened by ROC analysis. The MRL/lpr mice exhibited notable lymph node enlargement and characteristic rashes starting at 15 weeks of age (Fig. 5A). Additionally, an increased urinary albumin/creatinine ratio (UACR) was observed using the Bradford method. At 25 weeks, both groups were sacrificed, and plasma and kidney samples were collected for further analysis.. The results revealed that MRL/lpr mice showed significantly higher plasma ds-DNA concentrations and UACR compared to the C57BL/6 mice (Figs. 5B and 5C).

Through RT-qPCR detection, the five DEGs showed consistent changes with the analysis results of the GEO dataset. The mRNA expressions of PTPRC and MME were significantly different between the MRL/lpr group and the C57BL/6 group, and the expression of PTPRC was up-regulated while the expression of MME was down-regulated (Fig. 6A). Immunohistochemical staining analysis of renal histopathological sections also showed changes consistent with mRNA levels, and PTPRC, MME and MX1 also had significant differences at the protein level (Fig. 6B).