Identification of key biomarkers based on the proliferation of secondary hyperparathyroidism by bioinformatics analysis and machine learning

Human parathyroid tissue preparations

Parathyroid tissues were obtained between 2019 and 2022 from ten hemodialysis patients diagnosed with SHPT who underwent parathyroidectomy surgery, following approval from the Ethics Committee on Human Research of the Huashan Hospital, Fudan University, China (KY2018-369). Ultrasonography confirmed the presence of one or more PTGs >1 cm³. All participants signed an informed consent form after understanding their rights, the risks when participating in this study, and the purpose and process of our research. Immediately after resection, a piece of parathyroid tissue was taken from each gland, fixed with 4% paraformaldehyde, and embedded in paraffin for subsequent histological study. The rest tissues were frozen and stored at −80 °C. Each sample contains only one parathyroid tissue from one patient in subsequent immunohistochemistry and qRT-PCR testing. Basic demographics are shown in Table S1.

Animal model

All animal experiments were approved by the Ethics Committee for Animal Care and Use of Fudan University (No. 20171300A243). Male Sprague–Dawley rats (200–250 g) were purchased from Vital River Laboratory Animal Technology Co. Ltd. (Beijing, China). All the animals were housed in an environment with a temperature of 22 ± 1 °C, relative humidity of 50 ± 1%, and a light/dark cycle of 12/12 h. We randomly divided the animals into two groups of five animals each. A 5/6-nephrectomy (Nx) rat fed with high dietary phosphate is a commonly used animal model of human SHPT (Zhang et al., 2011). A 5/6-Nx was performed in a two-step procedure (Shvil et al., 1990). Following the second surgery, five rats were given a high-phosphate (1.2%) diet for 4 weeks. A control group consisted of five sham-operated rats fed with normal-phosphate diets (0.8%). Both the normal and high phosphate diets were obtained from Fanbo Biotechnology Co. (Wuxi, China) and contained the same levels of calcium (1.0%), vitamin D (750 IU/kg), lipids, and other nutrients. After a 4 weeks observation period, rats were anesthetized and sacrificed with sodium pentobarbital through intraperitoneal infection, and blood samples were obtained for creatinine, blood urea nitrogen (BUN), and intact PTH (iPTH). The PTGs were immediately removed by microdissection equipment. All rats were sampled in 1 day.

Data download

The GEO series matrix file GSE75886 was downloaded from the NCBI GEO public database, and the annotation file GPL10558. Thirteen patients of samples were included in the expression profile data; these consisted of five cases of hemodialysis patients with diffuse parathyroid hyperplasia, and eight cases of hemodialysis patients with nodular hyperplasia. We also downloaded the series matrix file GSE83421 and the annotation file GPL22020 for verification of key gene expression levels. In total, groups of samples were included in the expression profile analysis: six normal groups and 25 parathyroid adenoma patients.

Screening of key genes

The key genes identified from the GEO dataset were screened by differential analysis (p value < 0.05 and |logFC| >1), and the candidate gene sets were further screened by random forest (randomForest) and support vector machine (SVM) algorithms. Random forest is an ensemble learning algorithm with a decision tree as the base learner. It selects multiple samples from the sample set as a training set using the method of sampling and replacement, and uses the sample set obtained by sampling to generate a decision tree. At each generated node, it selects features randomly and without repetition, and then uses these features to divide the sample set, find the best dividing feature, and determine the predicted result. In this study, a random forest algorithm was used to evaluate the importance of candidate gene sets. A total of 1,000 classification trees were constructed, stirred 50 times, and evaluated according to the percentage of the increased mean-squared error (%IncMSE). SVM-Recursive Feature Elimination (SVM-RFE) is a machine learning method based on support vector regression analysis. SVM-RFE identifies the best variable (biomarker) by recursively removing variables to establish a more strongly supported classification model through the “e1071” software package which then further identifies the effect of these biomarkers on the disease. Following RF and SVM-RFE we kept the top-ranked features for subsequent analysis.

Functional annotation of differential gene sets

DEGs were annotated using the Metascape database (www.metascape.org). This allows a thorough exploration of the functional relevance of genes. Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genomes (KEGG) pathway analysis were also performed. A min overlaps greater than three and & pgreat were considered statistically significant. Min Overlap refers to the minimum allowable number or proportion of overlap between gene sets when comparing different gene lists. That is, if the number or proportion of overlap between two gene sets is lower than the set minimum threshold, it is not considered meaningful enrichment. Gene set enrichment analysis (GSEA) works based on a predefined set of genes. It sorts genes based on their differential expression between two groups, and then assesses whether the predefined gene set is more heavily associated with the top or bottom of the ranking list. In this study, GSEA was used to compare the differences in signaling pathways identified with KEGG and GO between the diffuse and nodular hyperplasia patients, and to explore the molecular mechanism of core genes between the two groups. The number of substitutions was set to 1,000, and the substitution type was set to phenotype.

Analysis of immune cell infiltration

The relative proportions of 22 types of immune infiltrating cells in patients in different subgroups were inferred through the analysis of gene expression data using the CIBERSORT algorithm. Pearson correlation analysis was performed on the risk score and immune cell content.

Regulatory network analysis of key genes

The R package “RcisTarget” was used to predict transcription factors. Calculations with RcisTarget are based on motifs, and the normalized enrichment score (NES) for these motifs depends on the total number annotated in the database. These motifs may be annotated either from the source annotated data file, or based on gene sequences similar to annotated motifs. To assess the overexpression of a motif in a dataset, the area under the curve (AUC) for each motif-motif set pair was calculated from the recovery curves of the gene sets for motif ordering. The NES for each motif was then calculated based on the AUC distribution of all motifs in the dataset. rcistarget.hg19.motifdb.cisbpont.500bp from the gene-motif rankings database was used.

Network prediction of lncRNA-miRNA-mRNA

Predict microRNA (miRNA) interaction pairs through the miRWalk database. In addition, in conjunction with TargetScan and miRDB databases, miRNAs shared by the three databases were selected as target miRNAs for key genes. Then, lncRNA was predicted through the ENCORI database, combined with lncRNA miRNA interactions and mRNA miRNA interactions to establish a lncRNA miRNA mRNA network, and visualized using Cytoscape.

Histologic studies

Part of the PTG tissues was embedded in paraffin and cut into 3 µm sections for hematoxylin and eosin (H&E) and immunohistochemistry (IHC) staining. The sections were stained using standard H&E and IHC protocols as previously described (Kan et al., 2018). According to the characteristics of HE staining, the parathyroid glands of SHPT patients were divided into diffuse and nodular hyperplasia (Mao et al., 2021). The primary antibodies used were CIDEA (1:100; Affinity, West Bridgford, UK), CAPZA1(1:50, Proteintech, Manchester, UK), and PCOLCE2 (1:100; Wuhan Fine Biotech, Hubei, China).

Quantitative real-time polymerase chain reaction (qRT-PCR)

Total RNA was isolated both from patient and rat parathyroid tissues using Trizol (Invitrogen) according to the manufacturer’s protocols and reverse transcribed with PrimeScript™ RT Master Mix (Takara). qRT-PCR was carried out with TB Green Premix Ex Taq (Takara). Primer sequences are given in Tables S2 and S3.

Laboratory measurements

Serum creatinine and blood urea nitrogen (BUN) levels were determined using an autoanalyzer (Siemens). Intact PTH (iPTH) was measured with rat-specific enzyme-linked immunosorbent assay (ELISA) kits (Immutopics).

Statistical analysis

All statistical analyses were performed in R studio (version 3.6). All results are expressed as mean ± SD. All statistical tests were two-sided and were performed using a t-test to determine p-values. Graphs were generated using GraphPad Prism 9. p < 0.05 was considered statistically significant.

Identification of key genes

Using the limma package to count the differential genes between the hyperparathyroidism nodular hyperplasia and diffuse hyperplasia samples in the GEO dataset GSE75886, with differential gene screening conditions p value < 0.05 and |logFC| >1, a total of 428 up-regulated genes and 186 down-regulated genes were screened out (Fig. 1A). These 614 differential genes were used as candidate genes for further analysis. Pathway analysis of candidate genes was performed through the Metascape database. The results showed that these candidate genes were enriched in lipid droplets, T-helper one cell differentiation, response to potassium ions, and other pathways (Fig. 1B). We then performed protein interaction network analysis on the genes in the candidate gene set with Cytoscape software (Fig. 1C).

The randomForest and SVM analysis further screened the 614 differential genes based on relation to disease disturbance (Figs. 2A and 2B). It was found that the sum of the ranking values obtained from the two analyses was distributed dispersedly (Fig. 2C). To more accurately screen genes related to the molecular mechanism of parathyroid glands, we performed a between-group rank-sum (Wilcox) test on the expression levels of the 614 differential genes in the validation set GSE83421; this dataset included six normal parathyroid glands and 25 adenomas from primary hyperparathyroidism (PHPT) patients. The key genes with a significant expression between groups and the best ranking in the machine learning process were identified as USP12, CIDEA, PCOLCE2, CAPZA1, and ACCN2 (Fig. 2D).

Analysis of key genes with immune infiltration

SHPT causes an increase in PTH release and the presence of PTH receptors on lymphocytes, leading to impaired B cell proliferation in patients with chronic kidney failure (Alexiewicz et al., 1990). There is an interactive relationship between SHPT and immune cells (Alexiewicz et al., 1990; Cozzolino et al., 2008; Komaba, Kakuta & Fukagawa, 2017), and when treating SHPT, it is necessary to pay attention to changes in immune cells to maintain the normal function of the immune system. The molecular mechanisms of the key genes affecting parathyroid progression were explored by looking at the relationship between these genes and immune infiltration in the parathyroid dataset. The proportion of immune cells in each patient was shown in Figs. 3A and 3B. When compared with the patients with diffuse parathyroid hyperplasia, the level of natural killer cells (NK) activated was significantly lower in the nodular hyperplasia group (Fig. 3C).

We further identified several key genes that were highly correlated with immune cells: CIDEA, PCOLCE2, and ACCN2 were all positively correlated with NK cell activation; USP12 and CAPZA1 showed a negative correlation with NK cell activation; and PCOLCE2 and ACCN2 were negatively correlated with macrophages M2 (Fig. 3D). We further observed correlations between these genes and different immune factors from the TISIDB database such as immunomodulators and chemokines (Fig. 3E). These analyses further support the close relationship between the identified key genes and immune cell infiltration levels and suggest that they may play important roles in the immune microenvironment.

GSEA analysis of key genes

The specific signaling pathways enriched by the five identified key genes and the potential molecular mechanisms of core genes affecting parathyroid progression were further explored with GSEA analysis. The USP12 gene GO-enriched pathways included “ncRNA export from nucleus” and “negative regulation of cell cycle arrest,” and the KEGG-enriched pathways included “endometrial cancer” and “glycosaminoglycan biosynthesis heparan sulfate,” among others (Fig. 4A). The CIDEA gene GO-enriched pathways included “interleukin-8 secretion” and “nuclear-transcribed mRNA catabolic process exonucleolytic,” and the KEGG-enriched pathways includes “aldosterone regulated sodium reabsorption” and “complement and coagulation cascades” (Fig. 4B). The PCOLCE2 gene GO-enriched pathways included “cell death in response to hydrogen peroxide” and “dendritic spine maintenance,” and the KEGG-enriched pathways included “aminoacyl tRNA biosynthesis” and “bladder cancer” (Fig. 4C). The GO-enriched pathways of the CAPZA1 gene included “de novo protein folding” and “positive regulation of chromatin organization,” and the KEGG-enriched pathways included “biosynthesis of unsaturated fatty acids” and “endometrial cancer” (Fig. 4D). The ACCN2 gene GO-enriched pathways included “definitive hemopoiesis” and “dendritic spine maintenance,” and the KEGG-enriched pathways included “biosynthesis of unsaturated fatty acids” and “drug metabolism other enzymes” (Fig. 4E).

Enriched motifs and corresponding transcription factors for key genes

The five identified key genes were found to be regulated by multiple common transcription factors which were identified by searching for enriched motifs using cumulative recovery curves and selection of important genes (Fig. 5A). The analysis showed that the motif with the highest normalized enrichment score (NES: 8.46) was annotated as cisbp_M6542, and the transcription factor ZBTB6 was identified as the master regulator. Two of the five genes were enriched in this motif: ACCN2 and USP12. All enriched motifs and corresponding transcription factors for the key genes are shown in Fig. 5B.

miRNAs and lncRNAs prediction of key genes

The prediction and analysis of the miRWalk database and ENCORI database were used to obtain the possible miRNAs and lncRNAs of these five key genes, respectively. First, the mRNA-related mRNA-miRNA relationship pairs were extracted from the miRWalk database, giving a total of 993 mRNA-miRNA relationship pairs of which 68 relationship pairs were retained (including three mRNA relationship pairs that were verified in TargetScan or miRDB and 68 miRNAs). Based on these 68 miRNAs, the interacting lncRNAs were predicted, resulting in a total of 4,152 predicted pairs of interactions (including 11 miRNAs and 2,224 lncRNAs). Finally, we constructed a ceRNA network in Cytoscape v3.7 (Fig. 5C).

Correlation between key gene and parathyroid-related disease-causing genes

A total of 2,333 parathyroid-related disease-causing genes were obtained from the GeneCards database (https://www.genecards.org/), and selected the top 20 genes with relevant scores to have a further analysis with the expression levels of the five key genes. The gene expression levels of key genes were found to be significantly correlated with the expression levels of multiple disease-related genes. Among them, CAPZA1 was positively correlated with CDKN1B (cor = 0.724), and USP12 was negatively correlated with GNAS (cor = −0.611) (Fig. 6A).

Predictive performance of key genes

The predictive performance of key genes was explored using the ROC curves of diagnostic efficacy validation; higher AUC values suggest better predictive performance. The AUC values of the five key genes were: USP12 = 0.8 (Fig. 6B), CIDEA = 0.85 (Fig. 6C), PCOLCE2 = 0.8 (Fig. 6D), CAPZA1 = 0.925 (Fig. 6E), and ACCN2 = 0.875 (Fig. 6F). These high AUC values suggest that these five key genes may be good predictors of the occurrence and development of disease.

Validation of key genes in different types of hyperplastic parathyroid tissues from uremic patients

To determine the expression of key genes in parathyroid glands in SHPT, qPCR was used to perform differential analysis between diffuse and nodular hyperplasia tissues. The results showed that gene expression of the five key genes (USP12, CIDEA, PCOLCE2, CAPZA1, and ACCN2) differed between diffuse and nodular hyperplasia tissues (Fig. 7A); CIDEA, PCOLCE2, and ACCN2 gene expression was significantly higher in diffuse hyperplasia tissues than in nodular hyperplasia samples. We further examined the expression of CIDEA, PCOLCE2, and CAPZA1 by IHC (Fig. 7B) and found that CAPZA1 was up-regulated in nodular hyperplasia of uremia patients, showing a more brownish patchy or granular staining when compared to that in the diffuse hyperplasia group.

Key Genes expression in animal models of SHPT

A 5/6-nephrectomy (Nx) rat fed with high dietary phosphate, as was done here, is a commonly used animal model of human SHPT (Zhang et al., 2011). Serum creatinine, BUN, and PTH levels were significantly higher in the Nx+HP group compared with the sham (Fig. 8A). The mRNA expression of CAPZA1 and USP12 in the parathyroid glands of SHPT rats was elevated, while that of CIDEA, PCOLCE2, and ACCN2 was decreased compared with sham rats. The results were similar to those seen in hyperplastic parathyroid tissues of SHPT patients (Fig. 8B).