A total of 8 PHPT patients (3 non-osteoporosis and 5 Osteoporosis) who underwent parathyroid resection surgery at Beijing Jishuitan Hospital were recruited in our study between January 2021 and February 2022. All patients underwent successful parathyroidectomy and were followed up from the time of diagnosis up to 36.0 months postoperatively. The diagnosis of PHPT was made mainly according to high or inappropriate PTH levels and the presence of hypercalcemia. Patients were included if they met the following criteria: (1) serum PTH level > 65 pg/mL and serum calcium level > 2.75 mmol/L; (2) parathyroid lesion excision performed by experienced physicians in the same department; (3) biochemical and BMD measurement before and after parathyroidectomy; and (4) patients diagnosed with symptomatic PHPT. Patients were excluded if they met the following criteria: (1) incomplete BMD measurements before and after parathyroidectomy or patients who could not be followed up; (2) normal parathyroid gland tissue (i.e. no hyperplasia, adenoma, and parathyroid cancer) diagnosed by histopathological examination after excision of the parathyroid lesions; and (3) serum calcium level remained above the normal range after excision of the parathyroid lesions. Signed informed consent forms were obtained from all subjects before the study. The resected parathyroid tissue samples were collected from all patients during the surgery and stored at -80°C. This study was reviewed and approved by the Institutional Review Board of Beijing Jishuitan Hospital (review batch number 201905–01).

In order to perform single cell sequencing data analysis, we followed both automated and manual procedures (15). For data loading and quality evaluation, we employed Seurat v4.0.1 package in R software (16). Following primary standards were adjusted to filter minimum level of cells; (i) total UMI counts < 1200, (ii) gene number < 300, and (iii) mitochondrial gene fraction > 20%. Based on aforementioned standards, in total 51624 cells including 24628 cells of non-osteoporosis and 26996 cells of osteoporosis patients suffering of osteoporosis were selected for further analysis. For data integration, Harmony package v0.1 was employed with default parameters (17). In total, 2500 high differentially expressed genes were identified and top 30 PCs were used for further dimensional reduction analysis. To ensure fidelity of parameters for cell clustering analysis, we determined resolution at 0.2.

For pseudo-time trajectory analysis in parathyroid cells, Monocle3 v1.2.9 and R package monocle v2.18.0 were individually employed (18), under default parameters. To discover root point, graph learning approach was employed.

In order to investigate cell communication and interaction, scRNA-seq data was analyzed with the help of CellChat v1.4.0 package in R software (19). In total 1,939 verified molecular interactions were considered in this study from CellChatDB (https://github.com/sqjin/CellChat).

Parathyroid chief cells secrete uncharacterized oxyphil cells in abundance in patients under treatment of hyperparathyroidism (20). We performed in-depth analysis to identify four different types of subpopulations of parathyroid cells in patients suffering with osteoporosis and non-osteoporosis and illustrated in UMAP. Sub-populations are comprised of two types of parathyroid chief cells (S100A13-PCC and CXCL10-PCC) and two types of oxyphilic cells (HSPA1A-OC and SPARCL-OC).

The functional enrichment analysis of Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) was conducted using ClusterProfiler v4.0 (21). We also evaluated the gene signatures scores using UCell v1.3 package (22), singscore v1.2.2 (23), AUCell v1.12.0 (24), GSVA v1.38.2 (25), and irGSEA v1.1.3 (26) in R software.

In order to reveal relative proportion of various cell types in parathyroid tissue, single cell sequencing data has been manually annotated into seven distinct cell types including parathyroid cells, fibroblast cells, T cells, endothelium cells, myeloid cells, mast cells, and B cells ( Figure 1A ). Among all, parathyroid cells displayed highest share 51.14% which shows their key role in osteoporosis as compared to all other type of cells, followed by 20.33% of endothelial cells, 11.19% of fibroblast cells, 7.93% of myeloid cells, 7.11% of T cells, 1.16% of B cells, and 1.14% of mast cells ( Figure 1B ). Similarly, proportion of parathyroid cells in non-osteoporosis patients was also higher as compared to patients suffering with osteoporosis ( Figures 1C, D ). Furthermore, classical markers expression analysis PTH revealed highest proportion of parathyroid cells among all types of cells, clearly depicted in Umap cluster ( Figure 1E ). Statistical analysis of relative proportion of different cell types in 5 patients suffering with osteoporosis and 3 non-osteoporosis patients were performed (27). We observed highest proportion of parathyroid cells in non-osteoporosis patients, while fibroblast cells and mast cells were significantly higher in patients suffering with osteoporosis (p < 0.05) ( Figure 1F ). Expression level of marker genes in each cell types is presented by feature plot and heat map ( Supplementary Figures S1A, B ).

We analyzed pathways and cell functions with the help of following four algorithms; AUCell, UCell, singscore, and ssgsea ( Figure 2A ). We observed no any significant variation in cell function and pathways between patients suffering with osteoporosis and non-osteoporosis. However, these algorisms displayed consistent variable trends on the base of cell types. Comparatively, myeloid, endothelial, and parathyroid cells displayed highest number of significantly up-regulated pathways ( Figure 2B ). These findings indicate that these four cell types are primarily involved in development of hyperthyroidism.

Pathway analysis shows that Hallmark-IL2-STAT5 and WNT-BETA-CATENIN signaling pathways were up-regulated, while oxidative phosphorylation pathway was down-regulated in osteoporosis patients as compared to non-osteoporosis patients ( Figures 2C, F, I ). Further investigation revealed that parathyroid, endothelial, fibroblast, and myeloid cells are involved in these variations ( Figures 2D, E, G, H, J, K ). These evidences further endorsed adaptability of these four types of cells, which probably play key role during pathophysiology of hyperthyroidism.

The cell function analysis was in agreement with cell communication analysis and vice versa. Among all cell types, highest interactions were observed among fibroblast, endothelial, parathyroid and myeloid cells ( Figure 3A ). Specifically, endothelial and fibroblast cells were in communication with parathyroid cells by secreting cytokines. The strength and number of cell interactions were higher in patients suffering with osteoporosis as compared to non-osteoporosis ( Figure 3B ). Moreover, differential interaction strength analysis revealed that fibroblast cells are hub cells playing key role in cell-to-cell interaction between both groups ( Figure 3C ). Further in osteoporosis patients, highest differential number of interactions and their strength was observed in fibroblast cells, which are heterogenic cell cluster and secrete highest signaling molecules to communicate with rest of the cells ( Figure 3D ). Fibroblast cells are predominantly involved in secretion of fibroblast growth factor (FGF) to affect parathyroid cells via FGF pathway ( Figures 3E, F ). Due to these reasons, patients suffering with osteoporosis differ in clinical symptoms as compared with non-osteoporosis.

Subcluster analysis of parathyroid cells revealed highest abundance 31.86% of CXCL10-PCC parathyroid chief cells, followed by 24.89% of HSPA1A-OC parathyroid oxyphil cells, 24.41% of SPARCL1-OC of parathyroid oxyphil cells and 18.84% of S100A13-PCC parathyroid chief cells in parathyroid tissue ( Figures 4A, B ). In patients suffering with osteoporosis, SPARCL1-OC followed by HSP1IA-OC parathyroid oxyphil cells were highly clustered, while in non-osteoporosis patients CXCL10-PCC followed by S100A13-PCC parathyroid chief cells were highly clustered ( Figures 4C, D ). GO enrichment analysis revealed that all four types of cells are functionally independent ( Figure 4E ).

Functional analysis revealed parathyroid cells in patients suffering with osteoporosis displayed highly up-regulated clusters as compared to non-osteoporosis ( Figure 5A ). We also found that the CXCL10-PCC, HSPA1A-OC, and SPARCL-OC considerably outperformed the S100A13-PCC in terms of up-regulated pathways ( Figure 5B ). In total, 281 differentially expressed genes (DEGs) were identified between non-osteoporosis and osteoporosis parathyroid cells including 177 up-regulated and 104 down-regulated genes ( Figure 5C ). GO analysis revealed that DEGs are highly enriched in regulation of kinase activity. PTH and CD59 were significantly upregulated in all kinds of parathyroid subpopulation in osteoporosis patients ( Figure 5B ). Contrarily, PTHLH was upregulated in non-osteoporosis patients.

Pseudotime trajectory analysis revealed that oxyphilic cells were more developed as compared to parathyroid chief cells ( Figure 6A ). In both non-osteoporosis and osteoporosis patients, cell communication analysis revealed that CXCL10-PCC, HSPA1A-OC, and SPARCL-OC are predominantly involved in cell-to-cell contact ( Figure 6B ). Compared with the non-osteoporosis, parathyroid cells in osteoporosis patients have significantly large number and higher strength of cell interaction ( Figure 6C ). Specifically, heat map analysis revealed that HSPA1A-OC parathyroid oxyphil cells play key role in large number and higher strength of interactions ( Figure 6D ). Different pathways interaction analysis revealed that osteoporosis patients had higher levels of NOTCH, PTH, FGF, and EGF pathway whereas non-osteoporosis patients had highest level of NRG pathway ( Figure 6E ). According to a PTH pathway interaction analysis, the SPARCL-OC subpopulation is the main source of PTH generation ( Figure 6F ).