Bioinformatics Analysis for Identification of Key Genes in Salivary Gland and the Potential of a Combination of Biomarkers for the Diagnosis of SS

Introduction

Sjogren’s syndrome (SS) is a systemic autoimmune disease primarily characterized by epithelitis in exocrine glands and the development of specific antibodies to SSA/SSB. Europe, North America, and Asia have a wide range of incidence rates for SS, and the prevalence rate in the population of three continents is 60.82 persons per 100,000 inhabitants,¹ and the female versus male ratio is 14:1. Patients with SS typically present with dry mouth and eyes, joint pain, and fatigue.^2,3 50% of patients present with systemic damage, and they have a 10–20-fold increased risk of developing B cell lymphoma.⁴

The etiology and pathogenesis of the characteristic progressive glandular dysfunction in SS have not been fully clarified. Genetic susceptibility, environmental factors, and immune disorders are thought to play a dominant role. Studies have found that susceptibility genes, such as human leukocyte antigen (HLA), Interferon Regulatory Factor 5 (IRF5), BLK, Interleukin 12A (IL12A), TNFAIP3 interaction protein 1 (TNIP1), and Chemokine receptor 5 (CXCR5), are closely related to the occurrence of SS.⁵ Activation of the type I interferon pathway, infiltration of T and B lymphocytes, and formation of ectopic germinal centers within salivary glands all cause damage in multiple tissues.⁶ The histopathology of salivary glands has the potential to both function as a biomarker and provide information regarding the prognostication and stratification in clinical trials.⁷

The main characteristic of SS upon biopsy, with 80.4% specificity and 86.2% sensitivity, is focal lymphocytic salivary adenitis (FLS).⁸ The key requirement of FLS is dense aggregates (foci) of ≥ 50 lymphocytes based on survey of at least four lobules, and the focus score (FS) calculation was defined as FS ≥ 1 focus per 4 mm² area.⁹ EULAR-SS task force clinical recommendations state that focal lymphocytic salialitis is considered as an important marker and can be used for the early diagnosis of SS.¹⁰ During the early stage of SS, there are no obvious signs of dryness, but the immune microenvironment of the salivary gland changes. Imaging biomarkers, such as ultra-high frequency ultrasound (UHFUS), may help reduce unnecessary and inadequate biopsies, but it has limited evidence of efficacy and is seldom used in the clinic.¹¹ With the rapid development of microarray technologies and bioinformatics analysis in the last few years, microarrays aimed at the genes and pathways of SS-like animals and SS patients have been studied.Nevertheless, single gene can not completely distinguish patients with from those with sicca symptoms or serological manifestations but no SS. Because of the highest item weights (scoring 3) of five items in 2016 American College of Rheumatology/European League Against Rheumatism classification criteria, variable sensitivity and specificity as well as low reproducibility, the diagnosis performance of MSGB could be further explored through multivariate and quantitative analyses of biomarkers.¹² Therefore, we analyzed three microarray datasets of salivary glands via bioinformatics analysis and further identified key genes and functions in SS salivary glands for higher sensitivity and specificity of MSGB, and we demonstrated the potential of multiplexed measurement of salivary protein biomarkers compared with FS using diagnosis model.

Materials and Methods

Microarray Datasets and the Identification of Differentially Expressed Genes (DEGs)

The salivary gland gene expression profiles of GSE23117 (n = 4 control, n = 11 SS), GSE7451 (n = 10 control, n = 10 SS), and GSE127952 (n = 6 control, n = 8 SS) were obtained from NCBI-GEO (https://www.ncbi.nlm.nih.gov/), which together included 20 healthy controls and 29 SS patients. We used GEO2R online software (http://www.ncbi.nlm.nih.gov/geo/geo2r) to analyze three microarray datasets, including GSE23117, GSE7451, and GSE127952, that were provided by the original submitter in order to determine the DEGs. A P value < 0.05 and log | FC | ≥ 1 were used as our standard screening criteria for DEGs. Genes with a log | FC | ≥ 1 were considered to be upregulated genes, and those with a log | FC | ≤ –1 were considered downregulated genes. We also used the online drawing tool (https://bioinfogp.cnb.csic.es/tools/venny/index.html) on the DEGs of the three datasets to construct Venn diagrams.

ClueGO and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis

ClueGO plugin in Cytoscape 3.7.2 was used to perform Gene Ontology (GO) functions, and we used the David database (https://david.ncifcrf.go) to obtain the KEGG enrichment analysis of the DEGs, and we visualized the obtained data through the R “Bioconductor” package. DEGs and signaling pathway analysis was performed with the use of GlueGO and KEGG.

Construction and Analysis of Protein–Protein Interaction (PPI) Networks

All upregulated genes were imported into the STRING database (http://www.string-db.org). The protein species was set to “HomoSapiens,” the lowest interaction threshold was set as “low confidence” (0.15), and a PPI network was imported into Cytoscape software for visual analysis. Network topology analysis was carried out to obtain the hub genes with higher degree value. Top 20 degree genes were calculated and identified by CytoHubba as SS preselected hub genes.

Validation of Hub Genes and Evaluation of Their Efficacy Using Hematoxylin and Eosin (H&E) and Immunohistochemistry

In order to confirm the expression of hub genes in SS patients and suspected patients who were not yet diagnosed, 28 labial gland biopsy samples, including SS patients and non-SS patients (Table 1), were obtained from patients admitted to The Second Affiliated Hospital, Zhejiang Chinese Medical University, Hangzhou, Zhejiang Province, P.R.C, who received a definitive diagnosis of SS (they were diagnosed by the 2016 American College of Rheumatology/European League Against Rheumatism classification criteria for primary Sjögren’s syndrome¹²). Informed consents were received from all enrolled patients prior to MSGB. The Second Affiliated Hospital, Zhejiang Chinese Medical University Ethics Committees approved the protocols, and the International Conference on Harmonization guidelines and the principles of the Declaration of Helsinki were followed. Labial gland samples were fixed, embedded, sectioned (4 μm), dewaxed with xylene, rehydrated through an alcohol gradient, and stained with H&E. A STAT1 Rabbit IgG Polyclonal antibody (10144-2-AP) from Proteintech (Wuhan Sanying, D3-3, No. 666, Gaoxin Avenue, Wuhan East Lake Hi-tech Development Zone, Wuhan, Hubei, P.R.C.) and MNDA Rabbit IgG Polyclonal antibody (DF4191), IL-10RA Rabbit IgG Polyclonal antibody (DF6643), and CCR1 Rabbit IgG Polyclonal antibody (DF2710) all from Affinity (Xiangtai Biological Technology Co., Ltd., No. 128, Shen Jiang Bian village, Changzhou, Jiangsu, P.R.C.) were used at a 1:200 dilution in this study. For immunohistochemical analysis, sections were processed for antigen retrieval with citrate buffer solution (pH 6.0) before blocking endogenous peroxidase activity with 3% H2O2 for 15 minutes. Sections were blocked with goat serum for 15 minutes and incubated for 1 hour at 37°C. To evaluate antigen expression, each section was photographed at ×200 magnification, and five images were taken across sections. The optical density was analyzed by Image-Pro Plus 6.0 (integrated optical density (IOD) = mean optical density × area). Receiver operating characteristic (ROC) curves were adopted to analyze the diagnostic value of the IOD of four genes for SS. The plot area under the curve (AUC) was computed through numerical integration of the ROC curves. The genes with the highest AUC value were identified as having the strongest power for diagnosing SS.

Table 1 Patients Clinical and Serological Information

Statistical Analysis

Bioinformatics statistical analysis was performed using R version 4.0.3, and GraphPad Prism version 8.0.2 (GraphPad Software, San Diego, CA, USA) was used for computing ROC analysis. Numerical results were presented as the means ± standard deviations (SD) or medians with interquartile ranges (IQR), and Wilcoxon Mann–Whitney U two-sample rank sum test was used to analyze the statistical significance between two groups. P<0.05 was considered statistically significant.

Results

Identification of DEGs in Salivary Glands

Three microarray datasets, which included SS patients and healthy controls, were analyzed with Venn tool. We identified 53 upregulated DEGs and 15 downregulated DEGs (data not shown) from the GSE23117, GSE7451, and GSE127952 datasets (Figure 1).

Figure 1 Venn diagrams show the DEGs in the salivary glands of patients with SS. Venn diagrams showing the upregulated DEGs (A) and downregulated DEGs (B).

ClueGO Enrichment Analysis

Enrichment analysis was performed on 53 intersecting targets using the ClueGO v2.5.7 plugin. After setting the P value as < 0.05, Kappa Score Threshold as 0.4, and Golevel as 3Minimum#Genes and 4%#Genes, we identified 31 GO terms. The results showed that the genes were mainly functionally enriched in dendritic cell migration, regulation of the extrinsic apoptotic signaling pathway via death domain receptors, lymphocyte homeostasis, regulation of B cell proliferation, and response to type I interferon (Figure 2).

Figure 2 ClueGO enrichment analysis. Significantly enriched GO terms of DEGs in SS based on their GlueGO (Biological Processes) and CluePedia functions (A). The pie chart of ClueGO enrichment analysis (B).

The Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Enrichment Analysis

The results of KEGG analysis showed that the upregulated genes were associated with several pathways, including the NOD-like receptor signaling pathway, human cytomegalovirus infection, cytokine–cytokine receptor interaction, and herpes simplex virus 1 infection (Figure 3A). Downregulated genes showed low correlation with SS because of low score compared with upregulated genes (Figure 3B).

Figure 3 KEGG pathway terms of DEGs in SS. A histogram of KEGG enrichment analysis of upregulated genes, which shows a gradual change in color from red to purple, indicating the change in the P value from small to large (A). A histogram of KEGG enrichment analysis of downregulated genes (B).

Construction of PPI Networks of the Upregulated Genes and Top 20 Hub Genes

PPI networks of the DEGs were constructed using STRING and Cytoscape. After removing three loose nodes, the PPI network contained a total of 50 nodes and 507 edges. We identified the top 20 hub genes, including PTPRC, CCR1, FGL2, MNDA, TNFSF13B, STAT1, IL10RA, LCP2, GBP5, GBP2, SAMD9L, NLRC5, TAGAP, CD53, IKZF1, GBP1, IFIT3, CXCR4, DOCK8, and XAF1 using Cytohubba according to degree. We found that the degree values were largest for Receptor-type tyrosine-protein phosphatase C (PTPRC), Signal transducer and activator of transcription 1(STAT1), Myeloid cell nuclear differentiation antigen (MNDA), Interleukin-10 receptor alpha (IL10RA), and C-C chemokine receptor type 1(CCR1) (Figure 4).

Figure 4 Protein–protein interaction (PPI) networks for 53 upregulated DEGs (A) and the top 20 hub genes (B) were constructed using Cytoscape visualization, which shows a gradual change in color from purple to yellow, indicating the change degree value from small to large.

Expression of Hub Genes in Affected Labial Gland Tissues

In order to validate the diagnostic value of the hub genes, H&E and immunohistochemistry analysis of four upregulated candidate genes, including STAT1, MNDA, IL10RA, and CCR1, in 14 patients and 14 controls were assessed. There was no obvious atrophy in the glands of non-SS comparator group, and the stroma was infiltrated with scattered lymphocytes and plasma cells. The glands of SS patients exhibited various degrees of acinar atrophy and destruction, ductal dilatation, and interstitial fibrosis as compared to non-SS patients, some SS glands showed granulomatous inflammation, adiposis, and germinal centers as compared to non-SS patients. STAT1, MNDA, IL10RA, and CCR1 proteins had increased expression in the labial gland of SS patients as compared to non-SS patients (Figure 5). The IOD of STAT1 and IL10RA was significantly higher in MSGB from SS patients than from non-SS patients (P < 0.001).

Figure 5 H&E and immunohistochemical analysis (A) of STAT1, MNDA, IL10RA, and CCR1 proteins in minor salivary gland biopsy (MSGB) staining from 14 SS patients and 14 non-SS patients (magnification: ×200 for H&E staining and ×200 for IHC). Blue marks the nucleus, and brown marks the target protein. The bar graph indicated medians with interquartile ranges of IOD value (B). ***P<0.001.

Combination ROC Analysis

The IOD obtained from 14 SS and 14 non-SS patients was converted by logistic regression, and then the ROC value of FS and these hub genes was calculated. Each AUC and combined predicted probability were computed. FS showed a good performance with an AUC of 0.893 (Figure 6A). STAT1 (AUC = 0.807; P < 0.001) and IL10RA (AUC = 0.818; P < 0.001) had a high performance, while CCR1 (AUC = 0.587; P < 0.001) had a low performance. The multiple-gene predicted probability (AUC > 0.9) was higher than the single-gene predicted probability and FS, which signified that the above four genes can be used to distinguish SS from non-SS as compared to MSGB (Figure 6B).

Figure 6 The single and combined diagnostic performance of the four genes. The diagnostic performance of FS (A). The diagnostic performance based on the IOD of four genes in SS patients (B). An AUC > 0.8 indicated that the predicted model had good efficacy. Abbreviations: ROC, receiver operating characteristic; AUC, area under the ROC curve.

Discussion

Although SS is one of the autoimmune diseases, the early diagnosis and treatment of SS are far less studied than systemic lupus erythematosus and rheumatoid arthritis.¹³ A better understanding of the pathological processes is still needed.¹⁴ MSGB classification criteria have been based on three main classification systems including Chisholm and Mason, which described the degree of T or B lymphocyte infiltration or focus of lymphocytes around the salivary glands, with a sensibility and specificity of 72.1% and 80%, respectively.¹⁵ Greenspan and Daniels, defined the concept of the FS and FLS, and Tarpley, which added concepts of acinar destruction and fibrosis. While MSGB facilitates a diagnosis in some seronegative patients, it only has a moderate diagnostic performance and needs to be further refined. Subtle changes in foci size might not be accurately described, and the sensitivity and specificity of these scoring systems have never been compared before.¹⁶ There are some patients who have atypical clinical presentations and an FS of zero. The further quantification of foci per 4 mm2 tissue section utilizing key genes expression will help to accurately stratify those patients whose clinical information is not obvious, but have a higher risk of developing as SS.¹⁷ Our study identified hub genes through the bioinformatics analysis of three salivary gland datasets from SS patients and healthy controls. The hub genes were validated via MSGB of SS patients and non-SS patients. The ROC of the combined diagnostic model using the four indicators had better performance compared with single indicator or FS.

Our ClueGO and KEGG analysis revealed that genes were mainly enriched in viral infection and interferon and inflammation-related signaling pathways, including the NOD-like receptor signaling pathway and human cytomegalovirus (HCMV) infection. Gland invasion and damage caused by viral infection have been implied as mechanisms of SS. HCMV establishes latency in early myeloid lineage cells, and interferon-responsive genes which are targeted during the establishment of HCMV latency. Patients with compromised immune function will relapse once HCMV reactivation occurs.¹⁸ A previous study found that SS patients had a significantly higher prevalence of IgG-anti-EA and IgG-anti-VCA antibodies (which are specific in population with acute EBV infections) than healthy controls, and patients with SS that had higher levels of IgG-anti-EA antibodies, which were linked to low levels of C3 and C4, in a case-control study.¹⁹ Interferon (IFN) plays an important role in both the pathophysiology and clinical manifestations of SS.²⁰ Most SS patients express a “type I IFN signature” in immune cells, which drive the production of auto-antibodies and a range of clinical manifestations.²¹ Infection and other factors are known to activate the expression of type I IFN and IFN-γ. These factors activate the NOD-like receptor signaling pathway or transcription of the STAT family, which induce expression of IFN-stimulated genes (ISGs), as well as feedback inhibitory molecules, such as IL-10, in diseased tissues.²² A previous study showed that increased virus-like repetitive elements, such L1 genomic repeat elements, represent endogenous nucleic acid triggers of the type I IFN pathway in SS and may contribute to injury of salivary gland tissue.²³

In order to find potential biomarkers for SS, the expression of STAT1 in labial gland biopsy samples was assessed. Phosphorylation of STAT1, a type of signal transducer and transcription activator, mediates cellular inflammation responses in both type I and II IFN signaling pathways.²⁴ Phosphorylated STAT1 induces complex of transcription factors 3(ISGF3) to translocate into the nucleus, where it binds to the IFN stimulation response element (ISRE) to activate the transcription of ISGs, leaving the cell in an antiviral state. Upon stimulation of IFNα2b, phosphorylation of STAT1 differed significantly between SS patients and healthy controls. Type I and II IFN-induced gene expression was positively correlated with the IFN-α2b-induced phosphorylation of STAT1 in B cells, which was associated with the production of auto-antibodies. Targeting STAT1 may contribute to the treatment of patients with SS.²⁵ A JAK-1 inhibitor was thought to inhibit IFN-induced transcription genes and B lymphocyte stimulator (BLyS) dominating B-cell differentiation, homeostasis in human salivary glands, as well as regulate IFNα- and IFNγ-induced pSTAT-1 in human salivary gland epithelial cells (SGECs).²⁶

Myeloid cell nuclear differentiation antigen (MDNA), considered to be an IFN-inducible gene, is a member of the family of the pyrin and HIN domain containing (PYHIN) proteins, and it plays an important role in innate immunity, especially in human granulocytes and monocytes and earlier stage cells in the myeloid lineage.²⁷ MNDA plays an emerging role in inflammation and apoptosis in SS.²⁸ It has been reported as a worthwhile biomarker in PBMCs, which can lead to highly sensitive and specific detection for SS (AUC = 0.84).²⁹

The IL10 receptor α (IL10 Rα) and IL10 receptor β (IL10 Rβ) are members of the IFN receptor (IFNR) family. The IL10 signal is transmitted through two homodimers and a heteroprotein tetramer receptor complex consisting of two IL10 Rα and two IL10 Rβ. STAT signaling is induced downstream by phosphorylation of the cytoplasmic tail of the IL10 Rα and IL10 Rβ by JAK1 and Tyk2, respectively.³⁰ IL10 suppresses proinflammatory responses and sometimes plays a protective role in SS. Possible linkage of IL10 polymorphisms to SS susceptibility has a critical role in the clinical symptoms and serological biomarkers.³¹ Functional mutations in genes, such as IL10 and IL10RA, have also been reported to be associated with immune dysregulation.³² IL10-producing regulatory B (Breg) cells have been shown to inhibit the follicular helper CD4 T cells (Tfh) cell response by promoting STAT5 phosphorylation and exerting an influence on the therapeutic effects in SS.³³ Immunohistochemistry of salivary gland adipose tissue demonstrated that SS patients candidate genes leading to inflammation, such as IL10, IL10RA, and IFN production.³³

C-C chemokine receptor 1 (CCR1) is one of the dominant regulators of leukocyte recruitment to resting sites in inflammatory disease. A previous study showed that CCR1-deficient mice exhibited blood and spleen defects, in addition to depletion of eosinophil numbers in the lungs.³⁴ Another study reported that CCR1 antagonists had long-term anti-inflammatory effects and were effective in blocking splenocyte migration in vitro. As such, targeting chemokine receptors is a promising treatment strategy for chronic inflammation diseases.³⁵ The role the CCR receptor family plays in leukocyte recruitment to salivary glands is unclear; therefore, the diagnosis efficacy of CCR1 is not important in this study.

Conclusion

Our results showed that both STAT1 and IL10RA contributed to the sensitive and specific of SS compared with MSGB. Detection of these biomarkers in salivary glands will provide a secondary application of MSGB, and demonstrate the potential of a combined measurement of salivary gland genes for the diagnosis of patients with SS. To our knowledge, this is the first study to develop a gene diagnosis model in SS patients using quantitative data of IHC. Further in vivo and in vitro experiments investigating the epigenetic regulation of interferon genes in lymphocytes are promising, and it is of great significance to the construction of the SS diagnosis model to synchronize interferon genes and inflammatory cytokines at the molecular level, which will help to realize precision medicine approaches of SS as soon as possible. However, it should be noted that this study only covered deficiency of salivary glands. Furthermore, evaluating the diagnostic and prognostic value of combined biomarkers of MSGB in a large cohort of patients with SS is necessary.