Integrated Bioinformatics and Validation Reveal Potential Biomarkers Associated With Progression of Primary Sjögren’s Syndrome

Background Primary Sjögren’s syndrome (pSS) is a chronic systemic autoimmune disease of the exocrine glands characterized by specific pathological features. Previous studies have pointed out that salivary glands from pSS patients express a unique profile of cytokines, adhesion molecules, and chemokines compared to those from healthy controls. However, there is limited evidence supporting the utility of individual markers for different stages of pSS. This study aimed to explore potential biomarkers associated with pSS disease progression and analyze the associations between key genes and immune cells. Methods We combined our own RNA sequencing data with pSS datasets from the NCBI Gene Expression Omnibus (GEO) database to identify differentially expressed genes (DEGs) via bioinformatics analysis. Salivary gland biopsies were collected from 14 pSS patients, 6 non-pSS patients, and 6 controls. Histochemical staining and transmission electron micrographs (TEM) were performed to macroscopically and microscopically characterize morphological features of labial salivary glands in different disease stages. Then, we performed quantitative PCR to validate hub genes. Finally, we analyzed correlations between selected hub genes and immune cells using the CIBERSORT algorithm. Results We identified twenty-eight DEGs that were upregulated in pSS patients compared to healthy controls. These were mainly involved in immune-related pathways and infection-related pathways. According to the morphological features of minor salivary glands, severe interlobular and periductal lymphocytic infiltrates, acinar atrophy and collagen in the interstitium, nuclear shrinkage, and microscopic organelle swelling were observed with pSS disease progression. Hub genes based on above twenty-eight DEGs, including MS4A1, CD19, TCL1A, CCL19, CXCL9, CD3G, and CD3D, were selected as potential biomarkers and verified by RT-PCR. Expression of these genes was correlated with T follicular helper cells, memory B cells and M1 macrophages. Conclusion Using transcriptome sequencing and bioinformatics analysis combined with our clinical data, we identified seven key genes that have potential value for evaluating pSS severity.


INTRODUCTION
Primary Sjögren's syndrome (pSS) is a chronic systemic autoimmune disease of the primary exocrine glands characterized by specific pathological features (1,2). pSS often features immune-mediated destruction of exocrine tissues, such as salivary and lacrimal glands, resulting in the loss of saliva and tear production (3). Immunological abnormalities such as antinuclear antibodies (ANAs), antibodies against Ro/Sjögren's syndrome-related antigen A (SSA) and La/Sjögren's syndromerelated antigen B (SSB), and hypergammaglobulinemia are often detected in pSS patients by laboratory tests (4,5). pSS is estimated to affect between 0.03 and 4.8% of individuals across different regions and global populations and occurs in women approximately 9 to 20 times more frequently than in men (6)(7)(8). However, the pathogenesis and underlying mechanism of pSS are still unclear due to the heterogeneity of clinical phenotypes and complex pathogenetic mechanisms. Therefore, identifying disease-associated biomarkers is necessary to understand the complex progression of pSS.
The hallmark characteristic histological feature of pSS is the infiltration T cells, B cells, plasma cells, macrophages, mast cells, and dendritic cells into the lacrimal and salivary glands (9). Salivary glands from pSS patients express a unique profile of cytokines, adhesion molecules, and chemokines compared to healthy controls (10). Nocturne et al. revealed that T lymphocytes are key factors in the immunopathogenesis of pSS and that the activation of B cells accelerates disease progression and drives some disease-specific phenotypes (11). Current studies have increased the understanding of autoimmune disease complexity. Although the innate immune system has been demonstrated to play a significant role in early stages of pSS, the exact role played by regulatory cells and their related genes in SS is far from understood. Labial minor salivary gland (LSG) biopsy is widely used in the diagnosis of pSS due to its simplicity, familiarity, and acceptance, and plays an irreplaceable role in the established American-European Consensus Group classification criteria and the proposed American College of Rheumatology/ European League Against Rheumatism (ACR/EULAR) criteria (12). Evidence suggests that investigation of LSGs could reveal potential biomarkers and additional prognostic, stratification, or mechanistic insights in pSS (13).
Previous Literature has reported that flow cytometry or immunohistochemistry can be used to detect certain immune cell subsets and genes related to pSS in peripheral blood mononuclear cells (PBMCs) and salivary glands (14). Over the past decade, genetic studies have been a useful and effective tool for the identification of novel pathogenesis-supporting pathways. Nocturne et al. reviewed the first genome-wide association study (GWAS) and identified three major pathogenic steps in pSS: innate immune system activation, B-cell activation, and T-cell activation (11). Using RNA sequencing data and four salivary glands microarray datasets, Luo et al. revealed the potential roles of ICOS in pSS (15). Recently, single-cell RNA sequencing has been applied to PBMCs from pSS patients, to identify the immune cell subsets and susceptibility genes related to the pathogenesis of pSS (16). However, little is known about genes associated with disease progression. Overall, delineating the molecular events and related genes involved in different stages of pSS represents an unmet clinical need, as this information could improve disease diagnosis and inform future therapeutic target selection.
Herein, we performed a comprehensive investigation involving transcriptomic analyses, salivary histology, and validation of hub genes in salivary gland samples from different stages of pSS. First, guided by integrated bioinformatics analysis, we identified differentially expressed genes (DEGs) from both our own RNAseq data (minor salivary gland) and pSS datasets (both major and minor salivary glands) from the NCBI Gene Expression Omnibus (GEO) database. Second, histochemical staining and transmission electron micrographs (TEM) were performed to macroscopically and microscopically characterize the morphological features of minor salivary glands at different disease stages. Furthermore, we performed quantitative PCR to validate hub genes. Finally, we analyzed the correlations between selected hub genes and immune cells in the autoimmune microenvironment. Our findings provide further insight into the mechanisms underlying the progression of pSS and suggest that the hub genes are potential diagnostic biomarkers and therapeutic targets.

Data Download and Preprocessing
The "GEOquery" package in R software was used to download the pSS datasets GSE40611 (17), 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), and GSE159574 (15) was downloaded from the corresponding supplementary file from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). The microarray datasets, GSE40611 and GSE127952, contained 17 pSS and 10 control parotid tissue samples and eight pSS and six control minor salivary gland samples, respectively. For these datasets, the normalize between arrays function in the R limma package was used to remove the inter-batch difference effects, and the maximum value method was used to remove the duplicate gene names corresponding to multiple probes. The high-throughput sequencing count datasets, GSE154926 and GSE159574, contained 43 pSS and 7 control minor salivary gland samples, and 16 pSS and 13 non-pSS salivary gland samples, respectively. For these datasets, the R Bioconductor package DESeq2 was used for data standardization and to obtain the standardized matrix file.

Bioinformatics Analysis
The DEGs from GSE40611, GSE127952, and GSE154926 were identified using the limma and DEseq2 packages (|logFC| >1, p<0.05) and using as training cohorts. The "ClueGO" plug-in with integration of Gene Ontology (GO) terms in Cytoscape was used to construct separate biological process networks from the list of upregulated genes (18). Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed using the cluster profiler package in R (19). P-value <0.05 was set as the cutoff criterion for significance. Interactive relationships and protein-protein interaction (PPI) networks of the overlapping DEGs were evaluated using the STRING database (http://string-db.org) (20). Cytoscape software (http://www.cytoscape.org/) was used to construct and visualize a biological network of key DEGs (21). A receiver operating characteristic (ROC) curve was used to analyze diagnostic values, including the area under the curve (AUC), sensitivity, and specificity obtained from the GSE159574 dataset. CIBERSORT (http://cibersort.stanford.edu) (22) was used to calculate immune cell proportion and estimate the scores of 22 human immune cell types in pSS samples from the GSE154926 cohort (p<0.05). Spearman correlation analysis was performed to assess correlations between hub genes and immune cells.

RNA Sequencing Data
We selected three cases each of pSS and non-pSS and sequenced each group of labial minor salivary gland tissue pools separately. Total RNA of each group was extracted by TRIzol Reagent (Invitrogen, USA) and quantified and qualified by Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA), NanoDrop (Thermo Fisher Scientific Inc.), and 1% agarose gel. One microgram total RNA with RIN > 6.5 was used for following library preparation. Nextgeneration sequencing library preparations were constructed according to the manufacturer's protocol. Then the prepared libraries were subsequently multiplexed and loaded on an Illumina HiSeq instrument (Illumina, CA, USA). Sequencing was carried out using a 2 × 150 bp paired-end (PE) configuration. Image analysis and base calling were conducted by the HiSeq Control Software (HCS), OLB, and GAPipeline-1.6 (Illumina) on the HiSeq instrument. To remove technical sequences, pass filter data of fastq format were processed by Cutadpt (V1.9.1) to be high-quality clean data. Differential expression analysis used the DESeq2 Bioconductor package. Padj of genes were set at <0.05 to detect differential expressed ones. The raw data were uploaded to the NCBI SRA database, under accession the number PRJNA722690.

Clinical Samples
Labial salivary glands from a total of 14 pSS patients with pSS, 6 non-pSS subjects (female, 40-65 years old), and 6 healthy control patients (mucous cysts or lower lip trauma) were collected from Ruijin Hospital, Shanghai Jiao Tong University School of Medicine (Shanghai, China). All pSS patients fulfilled either the 2016 ACR/ EULAR classification criteria (23) or the 2012 ACR classification criteria (24) for pSS. Non-pSS meet the same diagnostic criteria as patients presenting with xerostomia and xerophthalmia, but do not meet the classification criteria for pSS (25). The labial salivary gland biopsies were divided into four groups based on lymphocyte infiltration levels (25): healthy control, non-pSS, pSS1 (lymphatic infiltration foci = 1), and pSS2 (lymphatic infiltrating foci ≥ 2). Tissue samples for histochemical staining were fixed immediately following resection and were embedded in paraffin after 24 h. The samples for transmission electron micrographs (TEM) were fixed with 2.5% glutaraldehyde in a 0.1 M phosphate-buffered saline (PBS) for 2 h at 4°C. Other biopsy samples were flash frozen in liquid nitrogen and stored at −80°C until RNA extraction.

Histochemical Staining
Hematoxylin and eosin (HE) staining, Alcian blue (AB)-Sirius red staining, and Alcian blue-Periodic Acid-Schiff (PAS) stained sections from labial salivary glands of both non-SS and pSS subjects were examined under a light microscope (Leica, DMLB, Leica Microsystems Wetzlar, Germany). Inflammatory infiltration (grading and scoring), fibrosis (diffuse or partial), acinar atrophy with ductal epithelial changes, and amyloid deposits were assessed by independent observers. Tissue damage arising from both ducts and acini, lymphocyte infiltrates, fibrosis, etc., were also analyzed (26,27).

Transmission Electron Micrographs (TEM)
All non-SS and pSS labial salivary glands samples for transmission electron microscopy were subjected to 2.5% glutaraldehyde fixation, 1% osmium tetroxide postfixation, and ethanol gradient dehydration. Two changes of 100% propylene oxido (P.O.) for 10 minutes each and finally into the embedding resin. The samples were sectioned by diamond knife (LEICA EM UC7) to 70-90nm, followed by electron stained with lead citrate, and visualized on a Transmission electron microscope (HITACHI H-7650).

RNA Extraction and Quantitative RT-PCR
RNA was extracted using Takara RR420 and RR036A isolation kits according to the manufacturer's instructions. First-strand cDNA synthesis was performed using PrimeScript reverse transcriptase (Takara RR420). qPCR was carried out using the Takara RR036A kit. b-Actin was used as an internal reference. The relative expression of target genes was calculated using 2 −DDCt method. A P value <0.05 was considered significant. Primer sequences are summarized in Table S1. All PCR reactions were conducted in triplicate.

Immunohistochemical Staining of Biomarkers in Labial Salivary Glands
Based on the RT-PCR results, we selected two established proteins representing T cells and B cells, for validation. After deparaffinization, paraffin-embedded non-pSS, pSS1, and pSS2 labial salivary gland sections were incubated in antigen retrieval solution (Dako, Denmark), followed by the appropriate blocking steps. Following CD3 (IR503, Dako, Denmark) and CD19 (IR656, Dako, Denmark) primary antibody staining, slides were incubated with HRPconjugated secondary antibody and developed with DAB (Dako, Denmark). Fourteen pSS labial salivary glands and 6 non-pSS tissues were chosen randomly. Images were captured under a light microscope (Olympus Bx51, Japan).

Statistical Analysis
Statistical analysis was performed using unpaired Student's t-test or one-way ANOVA in GraphPad Prism software. Data are presented as means ± SD. A P value < 0.05 was considered statistically significant.

Identification of Common DEGs in pSS
The workflow of this study is shown in Figure 1. To investigate the common DEGs involved in pSS, we downloaded three public pSS datasets, including major (parotid) salivary gland (GSE40611) and minor salivary gland (GSE127952, GSE154926) datasets. As seen in Figure 2A and Table S2, a total of 47 upregulated and 2 downregulated common DEGs (≥3 instances of intersection) were identified (|logFC| >1, p<0.05). Figure 2B indicates whether these DEGs were found in the non-pSS versus pSS samples from each of three datasets. To systematically analyze the relationships between the common DEGs, we constructed a PPI network using the STRING database and visualized the data using Cytoscape 3.4.0 ( Figure 2C). In addition, a DEG network was constructed using DEGREE algorithm analysis within the Cytoscape software ( Figure 2D). Here, the hub nodes, including CD3G, CD3D, CD2, B2M, CD19, MS4A1, CXCL9, CXCL13, IFI6, GZMK, et al., were considered hub genes in the pSS DEG list. The identification of these DEGs laid the foundation for the subsequent identification of potential pSS biomarkers.

Function and Pathway Enrichment Analysis of Common Upregulated DEGs
To explore the potential biological processes associated with pSS, we performed pathway analysis on the common upregulated DEGs using the ClueGO plug-in within the Cytoscape suite (18). As seen in Figure 3A, we observed enrichment of immune-related signaling pathways like chemokine signaling pathway, T cell receptor signaling pathway, antigen processing and presentation, natural killer cell-mediated cytotoxicity, and Th1 and Th2 cell differentiation. Viruses that lead to inflammation, deregulation of epithelial cells, and autoimmune responses have been implicated as one of the pathogenic drivers of pSS (9). In line with previous studies, our findings demonstrated that common DEGs were also involved in infection-related pathways such as Epstein-Barr virus infection and viral protein interactions with cytokines and cytokine receptors. The pathway-related genes are displayed in Figures 3B, C. Using the DAVID database, the top five GO terms related biological processes among DEGs were primarily associated with response to virus, immune response-activating cell surface, receptor signaling pathway, T cell activation, antigen receptor-mediated signal pathway, and lymphocyte differentiation ( Figures 4A, B and Table S3). The terms related to molecular functions were mainly involved in cytokine receptor binding, antigen binding, and  Table S4). In addition, there is significant correlation in external side of plasma membrane, chromosome, centromeric region, and endocytic vesicle membrane in relation to cellular component ( Figures 4E-F and Table S5). These genes could be related to multiple biological pathways orchestrating pSS pathogenesis.

Identifying Candidate Biomarkers of pSS
RNA-Seq was carried out to further validate the results of the bioinformatics analyses above. As shown in Figures 5A, B, 28 DEGs were identified as common to the GSE127952, GSE154926, GSE40611, and our RNA-Seq data and therefore may be useful as candidate pSS biomarkers. To further verify the accuracy of these selected DEGs, another pSS dataset (GSE159574) was used as testing cohort to generate ROC curve. As shown in Figures 5C-E

Histological Observation of pSS Progression
As a limited amount of information is available regarding the link between macroscopic and microscopic features of the disease progression, we performed morphological observations  to explore changes to the acini and ducts during various disease stages. Clinical characteristics of the pSS patients are shown in Table 1. The labial salivary gland samples were divided according to their lymphocyte infiltration levels for labial salivary gland biopsy (25). Macroscopically, as seen in Figure 6A, we could clearly observe slight infiltration (Grade I); moderate infiltration or less than one focus (Grade II); one focus, mainly located in the periductal region (Grade III); and more than one focus with heavy interlobular and periductal lymphocytic infiltration (Grade IV). In addition, Alcian blue-Sirius red staining revealed abundant acinar tissue (blue) with clear cellular structures as well as evenly distributed collagen (red) in the acinar interstitium (Grade II) of non-pSS samples. pSS lesions of increasing severity, the serous and mucous acini appeared atrophied, and the blue and blue-purple stained areas were significantly reduced. In Grade IV samples, we noted severe acinar atrophy, increased interstitial collagen, and the appearance offat infiltration ( Figure 6B). Alcian blue-PAS staining showed that under normal circumstances, acidic mucous substances were blue, a glycogen and neutral mucous substances were red, and mixed mucous substances were blue-purple or purple-blue. As pSS disease severity increased, serous and mucous acini appeared atrophied, and blue and blue-violet stained areas were significantly reduced ( Figure 6C). Representative histochemical staining of the labial salivary glands revealed pSS stage-specific histopathology stages and suggested that related genes may also change significantly in severe pSS patients.

Scanning Electron Microscopy (SEM)
Ultrastructural alterations were detected in the labial salivary glands of pSS and non-pSS. In non-pSS (Grade II) sections, the acinar cell structure was normal, and the glandular epithelial cells in the acinus were myoepithelial cells. Cytoskeletal components and their unique myofilaments could be seen. By SEM we observed nuclear pyknosis of the duct epithelial cells and obvious swelling of organelles in pSS samples (Grade III) compared with non-pSS samples (Grade II) (Figure 7). In severe (grade IV) pSS, we observed disordered ductal structures, disappearance of some basal cells, obvious nuclear pyknosis of the duct epithelial cells, fat infiltration into the cytoplasm, and swelling of intracytoplasmic organelles (Figure 8).
To our knowledge, this is the first study to report the microscopic labial gland changes during pSS disease progression.

Identifying Potential Biomarkers Associated With pSS Procession
Although recent advances in high-throughput sequencing have provided the opportunity to identify disease-related genes in blood or tissue samples, relatively little is known about genes associated with disease progression. To explore potential  biomarkers associated with pSS stages and disease progression, we collected labial minor salivary gland biopsy samples from pSS patients and divided them into four groups according to lymphocyte infiltration. Subsequently, we performed RT-PCR to validate the expression levels of the abovementioned 28 DEGs. Of these, we selected seven hub genes ( Figures 9B, D) that were significantly higher in pSS2 samples (p<0.01) than in control, non-pSS, and pSS1 samples, including MS4A1, CD19, TCL1A, CCL19, CXCL9, CD3G, and CD3D (over 10-fold increase) (data for other genes not shown). In accordance with the results of previous reports and our own RT-PCR, we observed increased T cell (CD3 + ) and B cell (CD19 + ) infiltration in pSS2 (Grade IV) samples ( Figures 9A, C).

Correlation Between Potential Biomarkers and Immune Cells
The above results demonstrated that the hub genes were highly enriched in immune-related pathways and infection-related pathways. To explore the underlying mechanisms associated with these potential biomarkers, we extracted the gene expression matrix from the GSE154926 pSS samples and used the CIBERSORT algorithm to estimate correlations between these genes and infiltration of pSS samples with 22 human immune cell types. As shown in Figure 10, the potential biomarkers MS4A1, CD19, TCL1A, CCL19, CXCL9, CD3G, and CD3D were positively correlated with T follicular helper cells (Tfh cells), memory B cells, and M1 macrophage infiltration. This result indicated a strong correlation between the seven potential biomarkers and the immune response.

DISCUSSION
PSS is one of the more prevalent autoimmune diseases. The clinical presentation of pSS can differ greatly, ranging from common sicca symptoms to systemic symptoms (such as arthralgias, fatigue, neuropathy, or vasculitis) (28). PSS is a disease characterized by lymphocytic-mediated exocrine gland inflammation (29)(30)(31); however, the underlying genetic mechanisms driving the pathogenesis of pSS remain unclear. Lymphocytic infiltration of the exocrine glands, which leads to impaired lacrimal and salivary secretion, is the main characteristic of pSS (32). Although diagnostic criteria such as disordered saliva secretion have been widely used for pSS classification, their sensitivity still demands further improvement (33). The lack of standardized criteria used for diagnosis and classification of the disease leads to difficulty interpreting epidemiologic studies. In addition, a previous study reported that pSS patients typically do not receive a final diagnosis until as much as 7 years after the onset of symptoms (34). Therefore, the identification of potential biomarkers for pSS diagnosis is critical. Previous studies have extensively investigated molecular biomarkers of pSS over the last several decades; however, these studies were conducted based on only one kind of gland (3,8,(35)(36)(37). Considering the complexity of underlying pSS occurrence and development, it is essential for us to consider the expression of genes from all salivary glands in combination with clinical verification to identify more convincing biomarkers. In this study, we combined public datasets from major and minor salivary gland samples to screen for common DEGs as potential biomarkers. A total of 47 upregulated common DEGs were identified as likely to be involved in pSS development. The impact of DEGs is often to initiate fundamental molecular changes. To better appreciate the underlying mechanism of pSS, functional enrichment analysis was performed. The upregulated common DEGs were enriched in KEGG terms associated with immune-related pathways and infection-related pathways, including the T cell receptor and B cell receptor signaling pathways. GO terms implicated in pSS pathogenesis include immune system process (38), T cell activation (39)(40)(41), B cell activation (42,43), and response to virus (44,45). Previously reported pSS-related genes, including HLA-DRA (8,46), CCL19 (47,48), CXCL9 (49), and CD3D (50), were enriched in our study, suggesting that the related common DEGs are potential biomarkers for pSS. RNA-Seq was subsequently carried out to explore diagnostic biomarkers. The 28 common DEGs were thought to be directly involved in pSS development.
Investigating progression-associated gene expression profiles might enhance our understanding of the mechanisms underlying pSS progression. The public datasets we utilized for our bioinformatics analyses included transcriptomic profiles of both the major and minor salivary glands. To further screen for novel potential biomarkers associated with pSS progression  Grade III  27  1:80+  +  -++  +  +  40  1:320+  ++  --+  +  31  1:80+  -+  -+  +  69  1:160+  ---+  +  56  1:160+  --++  +  +  40 1:  and diagnosis, we collected labial salivary gland biopsy samples from 14 pSS and 6 non-pSS patients. We then investigated acinar and ductal morphologies at both the macroscopic and microscopic levels at various disease stages. Macroscopically, the occurrence of foci (an aggregate >50 lymphocytes), lymphocytic infiltration of the duct walls, and several other abnormalities (such as scattered lymphocyte infiltrates, interstitial sclerosis, and dilatation and regression of ducts and acini) are important features for the classification of pSS (51,51). Infiltration of lymphocytes into the salivary gland, which leads to destruction and subsequent fibrotic changes, is typically observed in pSS patients (52,53). Compared with non-pSS samples (Grade I, II), pSS samples clearly displayed one (Grade III) or more (Grade IV) foci. Microscopically, the development of pSS correlated with severe interlobular and periductal lymphocytic infiltration, increased acinar atrophy, and interstitial collagen and nuclear shrinkage and organelle swelling. To our knowledge, this is the first report using TEM to investigate the microscopic changes in pSS labial glands and the correlation of these changes with disease progression. T lymphocytes are considered to be key factors in the immunopathogenesis of pSS, and the activation of B cells accelerates the disease process and produces some disease-specific manifestations (11). In accordance with the enriched KEGG signaling pathways identified here, immunohisto-chemistry revealed that B cells and T cells are significantly increased around the ducts and acini of pSS samples compared with non-pSS samples. T lymphocytes (16,54) and B lymphocytes (55,56) were predominantly found around the ducts and within the lymphocytic infiltrates of pSS biopsy samples. Activated T cells contribute to disease pathogenesis by producing proinflammatory cytokines, which leads to the establishment of a positive feedback loop (57). Our observations demonstrate that pSS progression is associated with T/B cells, which is in accordance with functional enrichment analysis results (GO:0042110, GO:0050852, GO:0042113, and GO:0050853). Therefore, we selected common DEGs related to T cells (CXCL9, CD3D, CD3G) and B cells (MS4A1, CD19, TCL1A, CCL19) as potential biomarkers and verified these by qRT-PCR. Furthermore, we estimated positive correlations between potential biomarkers and Tfh cells, memory B cells, and M1 macrophages using the CIBERSORT algorithm. We suggest that heightened immune activation might broadly influence inflammatory infiltrates, which might subsequently influence the expression of DEGs.
MS4A1 (58) and CD19 (59) are B-lymphocyte-specific membrane protein-coding genes related to immunodeficiency. These genes play significant roles in the regulation of cellular calcium influx, development, differentiation, and activation of B lymphocytes. Mudd et al. indicated that MS4A1 is expressed in subsets of T cells and exhibits significantly higher expression in nivolumab-bound T cells from human patients undergoing anti-PD-1 immunotherapy (60), which is in accordance with the results of our functional enrichment analysis. MS4A1 has also been identified in a study published by Luo et al. (15). CD19+ B cells are often used for the identification of pSS and have been widely verified as a diagnostic marker (61,62). CCL19 maybe involved in more sever forms of autoimmune diseases, such as rheumatoid arthritis (63). Carubbi et al. indicated that CCL19 is involved in pSS pathogenesis (64), and Liu et al. identified CCL19 as a potential biomarker of pSS disease progression (47). CXCL9 has been shown to take part in T-cell trafficking, and Aota et al. found that enhanced production of CXCL9 may contribute to anti-inflammatory functions in pSS lesions (49). TCL1A is associated with precursor T-cell acute lymphoblastic leukemia and prolymphocytic leukemia (65). CD3G/CD3D forms the T-cell receptor-CD3 complex. Defects in CD3G are associated with T cell immunodeficiency (66), while CD3D is involved in T-cell development and signal transduction (67). Moreover, TCL1A is a novel biomarker that has never been reported in pSS. MS4A1, CD19, TCL1A, CCL19, CXCL9, CD3G, and CD3D might be potential diagnostic biomarkers associated with the progression of pSS, which remains to be verified based on a larger sample.

CONCLUSION
Our study describes characteristic pathological and ultrastructural changes in the salivary glands of severe pSS patients, including serous and mucous acinar atrophy, shrinkage of glandular epithelial cell nuclei, and damage to ductal structures and the basement membrane. Furthermore, using transcriptomic sequencing and bioinformatics analysis combined with our clinical data, we identified seven hub genes (MS4A1, CD19, TCL1A, CCL19, CXCL9, CD3G, and CD3D) that have potential value for the evaluation of pSS severity and that were strongly correlated with the immune response in the autoimmune microenvironment. Overall, understanding the mechanisms of pSS disease progression may pave the way for novel therapeutic interventions.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: https://www.ncbi.nlm.nih. gov/bioproject/PRJNA722690/.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethics Committee of Ruijin Hospital, School of Medicine, Shanghai Jiaotong University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
LJ and YG designed the overall research strategy. NL, LL, MW, YL, JY, YW, HX, DL, and XF performed the experiments. NL and LJ analyzed the data and wrote the manuscript. All authors contributed to the article and approved the submitted version.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021. 697157/full#supplementary-material FIGURE 10 | Associations between hub genes and immune cells as defined by the CIBERSORT algorithm. The horizontal axis represents immune cell contents, and the vertical axis represents genes. Different colors represent the correlation coefficient (red represents positive correlation, blue represents negative correlation), with a darker color indicating a stronger correlation. *p < 0.05, **p < 0.01, the asterisk represents the degree of importance.