LSG biopsies from 48 patients with pSS (without any related treatment) and 12 age- and gender-matched non-pSS sicca controls were performed in this study. The diagnosis of pSS was fulfilled according to the 2016 American College of Rheumatology/European League Against Rheumatism classification criteria (20) or the 2012 ACR classification criteria (21). Patients with non-pSS met the same diagnostic criteria as patients presenting with xerostomia and xerophthalmia but did not meet the classification criteria for pSS. Clinical information and samples were collected after patients signed written informed consent. The complete details are shown in Supplementary Table S1 . The Ethics Committee of Ruijin Hospital, Shanghai Jiao Tong University School of Medicine and Chinese Clinical Trial Registry, approved the study (ChiCTR2000039820). The patients with pSS were further stratified into two distinct stage groups, according to the severity of lymphocyte infiltration foci in labial salivary gland biopsy (13, 22). Patients with mild lesions (focal lymphocytic sialadenitis (FLS), with focus scores (FS) <2) were included in the low-infiltration stage (pSS1), and patients with severe lesions (FLS, with FS ≥2) were included in the high-infiltration stage (pSS2).

For histological staining, LSG samples were fixed freshly in 4% neutral formaldehyde overnight and embedded in paraffin. Samples were cut into 5-μm-thick serial sections. Immunohistochemical staining (IHC) was performed according to the manufacturer’s instructions. Briefly, slides were deparaffinized and microwave heated in citrate buffer (pH 6.0) to retrieve antigen. After gradual chilling, endogenous peroxidase activity was quenched using 3% hydrogen peroxide. Protein blockage was applied using 3% BSA for 30 min before incubation with primary antibodies at 4°C overnight. After washing with PBS, slides were incubated with secondary antibodies for 50 min at room temperature. Slides underwent color development with DAB (K5007, Dako, Denmark) followed by counterstaining in hematoxylin. The following primary antibodies were used: cytochrome c (1:500, ab133504, Abcam, UK) and Bcl-2 (IR614, Dako, Denmark). Finally, the slides were visualized under a light microscope (Nikon Eclipse Ni-U, Japan), and the images were captured using a camera attached to the microscope.

The ultrastructure of the minor salivary glands from non-pSS and pSS patients was visualized by transmission electron microscopy (TEM). Following fixation with 2.5% glutaraldehyde, the samples were postfixed in 1% osmium tetroxide and dehydrated using a gradient series of ethyl alcohol. Samples were then embedded in Embed 812 resin (EMS, TED PELLA, USA) and propylene oxide solutions followed by embedding in embedding resin for 48 h. The blocks were sectioned transversely at 70–90 nm using a diamond knife (EM UC7; Leica, Wetzlar, Germany). Ultrathin sections were stained with lead citrate and photographed with a transmission electron microscope (H-7650; Hitachi, Tokyo, Japan).

The integration of our clinical data and bioinformatic analyses is illustrated by the flowchart in Figure 1 . The pSS cohorts with publicly available datasets were obtained from GEO databases (23), including GSE40611 (24), GSE127952 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE127952), and GSE154926 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE154926). For validation, RNA-seq data and clinical information of an additional 39 patients with pSS were obtained from another publicly available dataset GSE173808 (25). We stratified the patients into three distinct groups: non-pSS (n = 12 labial glands, n = 14 parotid glands), pSS-low infiltration (focus score (FS) <2, n = 18 labial glands, n = 14 parotid glands), and pSS-high infiltration (FS ≥2, n = 14 labial glands, n = 7 parotid glands) (13). The Limma package and “DESeq2” package of R v3.6.1 were used for array data sets and high-throughput sequencing count data standardization, respectively, and the standardized matrix file is obtained. When multiple transcript IDs were present for one gene, we chose the ID with the highest average expression. Raw gene expression data were log2 transformed and quantile-normalized for all subsequent downstream analyses.

DEGs were selected using the “limma” and “DEseq2” R packages with a maximum posteriori absolute log2|fold-change| ≥1 and a Benjamini–Hochberg p-value <0.05. The Integrated Mitochondrial Protein Index (IMPI) of the MitoMiner database (http://mitominer.mrc-mbu.cam.ac.uk/) provides 1626 human genes that encode mitochondrially localized proteins (26). Overlapping mDEGs based on the MitoMiner database were extracted from GSE40611, GSE127952, and GSE154926, respectively, and visualized with a heatmap, Venn diagram, and volcano plot. The list of DEGs was used for GO (http://geneontology.org/docs/go-citation-policy/) and KEGG enrichment analyses (https://www.genome.jp/kegg/kegg1.html) (27), using the clusterprofiler package of R software (28, 29). Bioinformatic pathway analysis was conducted with the Gene Set Enrichment Analysis (GSEA) (https://www.broadinstitute.org/gsea/). GSEA is a computational method to detect statistical significance, and pathways using the KEGG gene set (c2.cp.kegg.v7.4.symbols.gmt) from the Molecular Signatures Database (MSigDB) (http://software.broadinstitute.org/gsea/msigdb/) by the JAVA program were selected as the reference gene sets (30). The algorithm of random sampling was 1,000 permutations. Only enrichment pathways with p < 0.05 and false discovery rate < 0.05 were considered statistically significant.

Heatmaps of the mDEGs were generated by “pheatmap” package v1.0.8 (https://CRAN.R-project.org/package=pheatmap) of R. To further investigate the correlation between hub genes and immune cell infiltration, formatted data were uploaded to the Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) R program (http://cibersort.stanford.edu). The CIBERSORT algorithm is a deconvolution algorithm that has been validated on gene expression profiles measured by RNA-sequencing, and it derives a p-value for the deconvolution of each sample with p < 0·05 considered accurate (31). We used CIBERSORT to estimate the fraction of 22 immune cell types in salivary glands from patients with pSS (GSE154926). The correlation between each hub gene and the 22 immune cells was tested by Spearman’s rank correlation and presented as a lollipop chart. We also performed functional enrichment analyses using the KEGG pathway dataset from GSEA.

Marker genes for immune cell types were identified (32). Of 24 gene signatures, 11 [e.g., dendritic cells (DCs), eosinophils, mast cells, macrophages, neutrophils, and natural killer cells (NKs)] were for immune cells in adaptive immunity, 13 (e.g., B, CD8⁺ T, T helper 1 [T_h1], T_h2, central memory T [T_cm], effector memory T [T_em], regulatory T [T_reg], and T follicular helper [T_fh] cells) for innate immunity. The immune infiltration score was calculated using the single-sample GSEA (ssGSEA) method (33) based on the GSE173808 dataset. Changes in the immune score between high and low hub gene expression subgroups were analyzed using the “ggpubr” package (https://github.com/kassambara/ggpubr) via a Wilcox test. Heatmaps and clustering analyses were generated using the “ComplexHeatmap” v2.10.0 package in R to show the correlation (34). Spearman correlation based on the genes and immune cells of interest was performed using the “pheatmap” package.

From the MitoCarta3.0 database (35) and high-throughput sequencing data (GSE173808), the DEGs of the oxidative respiratory chain complex were computed using the “limma” package and visualized with the “pheatmap” R package. The correlation between the four hub genes and five oxidative respiratory chain complex genes was calculated using Spearman’s rank correlation and visualized using the “ggplot2” R package (https://cran.r-project.org/web/packages/ggplot2/ggplot2.pdf) (36). We assessed the interrelationship between four hub genes and selected OXPHOS genes using Pearson’s correlation (R) with the function scatter from the “ggpubr” R package v0.4.0 (https://CRAN.R-project.org/package=ggpubr).

Correlations between four hub genes and mitochondrial metabolism, damage-associated molecular patterns (DAMPs) were computed with the Mantel test (37) and the Pearson correlation coefficient in pSS-low-infiltration and pSS-high-infiltration groups. The “ggcor” R package v0.9.8.1 (https://github.com/houyunhuang/ggcor), based on “ggplot2,” was used to provide a graphical display of any correlations and their combinations. In addition, we analyzed the expression of DAMP-related genes (NLRP3, ZBP1, TNF) and apoptosis-related genes (Bcl-2, Bax, caspase3) using the “ggpubr” package via a Wilcox test.

The LSG samples collected from patients were immediately immersed in the Allprotect™ Nucleic Acid and Protein Stabilization Reagent (R0121, Beyotime, Shanghai, China). Then, RNA was extracted and first-strand cDNA synthesis was performed using PrimeScript™ RT Master Mix (No. RR036A, Takara, Shiga, Japan), and qPCR was performed by TB Green^® Premix Ex Taq™ II (No. RR420A, Takara, Shiga, Japan). Primer sequences are summarized in Supplemental Table S2 , and β-actin was applied as an internal reference. The relative expression of target genes was calculated by the 2^−ΔΔCt method. All PCR reactions were conducted in triplicate.

Quantitative result data are presented as means ± SD. Statistical analysis was carried out with Student’s t-test or one-way analysis of variance in GraphPad Prism software. Statistical significance was set at p <0.05.

In this study, three publicly available datasets—GSE40611, GSE127952, and GSE154926—which contained 17, 8, and 43 patients with pSS were used as training datasets. In addition, the GSE173808 dataset which includes 39 patients with pSS was applied as a validation cohort. We then analyzed the DEGs associated with mitochondria based on the MitoMiner database. Heatmaps representing the most significant mDEGs (log2|fold-change| >1) are shown in Figure 2A and revealed a clear distinction between patients with or without pSS. Figures 2B, C show Venn diagrams representing the overlap between these mDEGs. A volcano plot represents the mDEGs (≥2 times intersection) between patients with pSS and controls. The top upregulated DEGs included IFIT3, CMPK2, and CD38 ( Figure 2D ).

Functional analysis of upregulated mDEGs revealed enrichment in KEGG pathways related to metabolism, translation, cell growth and death, ribosome, thermogenesis, and OXPHOS in patients with pSS ( Figure 2E ). The downregulated KEGG pathways involved protein processing in the endoplasmic reticulum, biosynthesis of amino acids, carbon metabolism, and metabolic pathways ( Figure 2F ). In addition, three domains of gene ontology (GO; biological process, molecular function, and cellular component) were analyzed using the GO database. The mDEGs upregulated for pSS were related to the biological processes OXPHOS, nucleoside triphosphate metabolic process, and adenosine triphosphate (ATP) metabolic process. Among the most relevant downregulated BP terms were that of the small-molecule catabolic process, co-enzyme metabolic process, electron transport chain, and fatty acid beta-oxidation ( Figures 2G, H ). The complete details are shown in Supplementary Table S3 .

The patients from the GSE154926 dataset were further used to screen out a mitochondrial-related gene signature. The expression heatmaps were generated by the “pheatmap” R package, and 21 mitochondrial-related genes (CD38, CMPK2, ITIF3, LAP3, TBC1D9, XAF1, IFI6, PMAIP1, IFI27, COASY, DNAJC4, GPT2, ITGA3, PYCR1, SERHL2, NME4, NT5DC2, SLC25A29, OGDHL, P4HB, SPSB3) were identified (p < 0.05). Among them, eight genes (p < 0.001) were further verified by real-time PCR ( Figure 3A ). The results suggested that the genes CD38, CMPK2, and TBC1D9 were relatively overexpressed in LSGs from all patients with pSS, while PYCR1 was underexpressed (p < 0.05). To further investigate the relationship between hub genes and infiltration of immune cells, the CIBERSORT algorithm was used on data of patients with pSS. As shown in Figures 3B, C , PYCR1 had a significant negative correlation with infiltration of memory B cells, M1 macrophages, CD8⁺ T cells, and T_fh (p < 0.05), while CD38 and CMPK2 had an opposite trend (p < 0.05). In addition, TBC1D9 had a positive correlation with infiltration of resting DCs, M1 macrophages, and T_fh but a negative correlation with plasma cells (p <0.05). In contrast, PYCR1 showed an opposite trend (p < 0.05). Subsequently, a total of four mitochondria-related genes (CD38, CMPK2, TBC1D9, PYCR1) were identified as the most promising factors associated with pSS disease severity.

GSEA was used to obtain a deeper insight into the function of DEGs. As shown in Figure 3D , overexpression of hub DEGs (CD38, CDMK2, TBC19) was highly enriched in pathways related to T cell/B cell receptor signaling pathways, primary immunodeficiency, NK cell-mediated cytotoxicity, cytokine–cytokine receptor interaction, and the JAK/STAT signaling pathway (p < 0.05). The pathways altered by PYCR1 were involved in the T cell receptor signaling pathway (p < 0.05, FDR 0.079), cytokine–cytokine receptor interaction (p < 0.01, FDR 0.049), and the JAK/STAT signaling pathway (p < 0.01, FDR 0.051) ( Supplemental Table S4 ). Interestingly, the GSEA results confirmed a strong association between the above four mitochondria-related genes and immune-related signaling pathways.

The relative level of immune cell infiltration for each patient (data from publicly available datasets mentioned above) was investigated for pSS. We compared the signature score of 24 immune cells using ssGSEA in patients stratified by histological focus score (non-pSS (FS = 0) vs. pSS-low infiltration (FS < 2) vs. pSS-high infiltration (FS ≥ 2)) depending on the validation cohort. The Kruskal–Wallis test revealed that the DC family (activated DC, immature DC, DC, plasmacytoid DC [pDC]), B cells, CD8⁺ T cells, cytotoxic cells, NK cells, T_h1/2 cells, T_fh, and T_reg cells were significantly higher in the pSS-high-infiltration group (p < 0.05) ( Figure 4A ). The bar plots in Figure 4B show the proportion of 22 immune cells in the pSS1 and pSS2 groups. The six most common immune cells in the pSS-high filtration group were memory B cells (26.7%), resting memory T cells CD4⁺ (15.2%), plasma cells (13.7%), M2 macrophages (10.3%), CD8⁺ T cells (9.2%), and M0 macrophages (7.8%), while the six most common immune cells in the pSS-low-filtration group were plasma cells (23.6%), resting memory B cells (17.4%), resting memory CD4⁺ T cells (16.0%), M2 macrophages (9.6%), naïve B cells (7.6%), and CD8⁺ T cells (7.0%). As shown in Figures 4C–F , there was some evidence of higher B cells, CD8⁺ cells, cytotoxic cells, T cells, T_cm, T_em, and T_reg immune scores for high CD38, CMPK2, and TBC1D9 expression (p < 0.001), while PYCR1 showed an opposite trend (p < 0.05, Wilcox test). In the high-immune cell infiltration group, CD38 was also positively correlated with B cells, CD8⁺ T cells, and cytotoxic cells, while PYCR1 was negatively correlated with cytotoxic cells, pDC, neutrophils, and T cells ( Figure 4G ).