The Institutional Review Board of the Gangnam Severance Hospital (IRB No. 3-2020-0160) approved this study. The sicca cohort comprised patients who visited the clinics of ENT, Eye, or Rheumatology in the Gangnam Severance Hospital with sicca symptoms, including xerostomia and xerophthalmia. The patients were administered a questionnaire assessing the following: “I have trouble swallowing solid food,” “I wake up while sleeping for a cup of water,” “I feel dryness in my mouth while eating,” “I need a sip of water to swallow solid food,” “I usually drink water to reduce dryness,” “I have trouble eating dry food,” “I feel dryness in my mouth while speaking,” and “I feel dryness in my mouth while chewing” (16). The patients with definite subjective dryness in their mouth or eyes were then enrolled in our sicca patient cohort. The saliva samples were collected for microbiome analysis from 43 patients. Of these 43 samples, 32 quality-controlled samples in which the concentration of the complementary DNA (cDNA) library was > 5 ng/μL, were used for the final analysis. Among the 32 samples, 23 were from patients with primary SS sicca who were diagnosed according to the recent ACR/EULAR classification criteria revised in 2016 (17), while nine were from non-SS sicca subjects who did not fulfill the diagnostic criteria. Patients with SS sicca exhibited positive anti-Ro/SSA results or a focus score in the LSGB sample of more than 1 foci/4 mm² with at least one of the aberrant items related to salivary or lacrimal glandular hypofunctions. Pathologists and radiologists who were blinded to the clinical and laboratory findings in this study analyzed the focus score in the LSGB samples and magnetic resonance (MR) sialography features of PGs. The exclusion criteria were as follows: history of irradiation of the head and neck area; radioactive iodine treatment for thyroid cancers; salivary gland surgery; other autoimmune diseases (potential secondary SS). Although some patients did not complete the radiological or laboratory examinations due to the coronavirus disease-2019 (COVID) pandemic, their clinical data were included for analysis once the statistical comparison was allowed.

The saliva samples were obtained from both PGs through salivary gland lavage. Patients were instructed to refrain from eating, drinking, and smoking for 2 h before salivary gland lavage. To reduce the chances of potential contamination, sample collection methods were performed as follows: The subjects rinsed their mouth with bottled water before sample collection. Citric acid (sialogogue) was administered to the floor of the mouth to stimulate salivation. After a small amount of saliva was secreted by PG massage, a rubber angiocatheter (Becton Dickinson Medicals Ltd, Singapore) connected to a syringe filled with 1 mL of physiological saline was gently inserted into the orifice of the PG at the opposite side of the upper second molar tooth. Gentle lavage was performed and saliva was regurgitated through the negative pressure of the syringe with PG massage. At least 0.1 mL saliva sample was collected from each side of the PG. Two syringes containing saliva samples from both sides were transported to the laboratory on dry ice. The saliva samples were aliquoted evenly into 1.5 mL Eppendorf tubes (Eppendorf, Hamburg, Germany) and stored at −80°C.

MR examinations were performed using a 3.0T MRI unit (Discovery 750; GE Healthcare, Milwaukee, WI) equipped with a quadrature head coil. The magnetic resonance imaging protocol was as follows. Axial T2 weighted IDEAL water imaging was performed using the following parameters: repetition time (TR), 7600 ms; effective echo time (TE), 85 ms; matrix, 320 × 256; field of view, 20 cm; section thickness, 3 mm; flip angle, 111; asymmetric echo shifts, −π/6, π/2, 7π/6. Axial T1 weighted images (T1WI) were obtained with the following parameters: TR, 780 ms; TE, minimal; matrix, 320 × 256; field of view, 20 cm; section thickness, 3 mm. MR sialography was performed using a three-dimensional fast-recovery fast spin-echo sequence in the axial plane (TR, 2000 ms; TE, 695.5 ms; echo train length, 110; field of view, 20 cm; matrix, 320 × 320, thickness, 1 mm). Salivation was stimulated by intraorally applying a sialogogue to enhance the visualization of the ductal structures. Sialography images were generated by creating a maximum-intensity projection. The data volume was reduced by eliminating the post peripheral sections with a high signal from the skin adjacent to the coil and removing the volume that included the eye and cerebrospinal fluid. Based on the MR sialography, structural deformation was assessed ( Figure 1 ). Sialectasis was defined based on the following Tonami’s criteria (18): Stage 1, punctuate with a diameter of ≤ 1 mm; Stage 2, globular with a diameter of 1–2 mm; Stage 3, cavity with a diameter of > 2 mm. Fat stage was evaluated based on the characteristic appearance of fat signals on T1WI and T2W IDEAL water imaging and using the grading methods reported by Izumi and Regier (19, 20). An iterative least-squares decomposition algorithm was employed to output a fat fraction map and an R2* map (21). A circular ROI was drawn within the PG. The ROIs were placed in the same position on the fat fraction map and R2* map.

Salivary flow tests and salivary gland scintigraphy were performed to determine the salivary gland function. A salivary flow test was conducted for 5 min by a nurse practitioner. The weight of the paper cup before and after the saliva was spit into it was measured. All patients were instructed not to eat, drink, or smoke for 2 h before saliva collection. For stimulated saliva collection, patients were administered 500 mg of ascorbic acid (sialogogue) for 4 min and instructed to continue spitting saliva for another 5 min into a second pre-weighed paper cup. The weight of the saliva was divided by 5 min to measure the salivary flow rates (SFRs) under unstimulated (UWMSFR, unstimulated whole mouth salivary flow rate) and stimulated (SWMSFR, stimulated whole mouth salivary flow rate) conditions.

Salivary gland scintigraphy was conducted using pertechnetate and Siemens gamma camera equipment as reported previously (22). The images were obtained for 20 min immediately after pertechnetate IV injection. Additionally, 5-min images were captured after stimulation with an oral sialogogue. At 1 min post-injection, the background count was checked until it reached the maximum count. The background count was the minimum when salivary secretion was stimulated after 20 min. The uptake ratio (UR) was calculated as the ratio of the maximum count to the background count. Tmin was defined as the time interval between the maximum and minimum counts. The maximum accumulation (MA) and maximum secretion (MS) were calculated as follows: MA = (maximum count − background count)/maximum count × 100 (%); MS = (maximum count − minimum count)/maximum count × 100 (%).

Saliva samples were stored at −80°C before sequencing. DNA was extracted from saliva samples using the DNeasy PowerSoil kit (Qiagen, Hilden, Germany), following the manufacturer’s instructions. The extracted DNA was quantified using Quant-IT PicoGreen (Invitrogen, Waltham, MA) and used for library construction according to the Illumina 16S metagenomic sequencing library protocols. To amplify the V3 and V4 regions of 16S rRNA, the following primers with Illumina adapter overhang sequences were used: V3-F: 5′-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG CCT ACG GGN GGC WGC AG-3’, V4-R: 5’-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA GGA CTA CHV GGG TAT CTA ATC C-3’. The purified PCR products were quantified using quantitative real-time polymerase chain reaction (qRT-PCR) according to the qPCR Quantification Protocol Guide (KAPA Library Quantification kits for Illumina Sequencing platforms). The amplicons were subjected to quality control using the TapeStation D1000 ScreenTape (Agilent Technologies, Waldbronn, Germany). Paired-end (300 bp each) sequencing was performed using the MiSeq™ platform (Illumina, San Diego, CA). All the above procedures were performed by Macrogen Inc. (Seoul, Republic of Korea).

ASVs in the paired-end sequence files were identified using the DADA2 pipeline (v1.21) (23). Briefly, sequence reads were filtered with expected error thresholds of 2 (start) and 5 (end) of each read. Filtered reads were de-replicated and denoised using the default settings in DADA2. After merging the paired reads and removing the chimeras, taxonomy was assigned using the Silva reference database (v138.1) (24). The Ape package (v5.5) was used to construct the phylogenetic tree. All the above analyses were performed in the R environment (v4.1).

Phyloseq objects were prepared with phyloseq (v1.37) and used for diversity analysis. Taxa that were only observed once in the entire dataset or appeared in one sample were excluded from further analysis. The observed ASVs, Shannon index, and Simpson’s diversity index were used to assess alpha diversity. Statistical analysis and data visualization were performed using Prism 7 software (GraphPad, San Diego, CA, USA). To compare the two groups, parametric t-test with Welch’s correction or nonparametric Mann-Whitney test was performed based on Gaussian distribution. Differences were considered significant at p < 0.05. Parametric one-way analysis of variance (ANOVA) or nonparametric Kruskal-Wallis test was performed based on Gaussian distribution for comparing more than three groups. The false discovery rate was calculated using the Benjamini-Hochberg procedure. An adjusted p-value of less than 0.05 was considered significant. Normality tests were performed using the D’Agostino-Pearson omnibus normality test or the Shapiro-Wilk normality test. The Bray-Curtis distance method was applied to assess beta diversity, and the results were plotted in ordination with principal coordinate analysis. A permutational ANOVA method was used to test the significance between groups based on the distance matrix. All statistical analyses were performed using the ordinate function in the phyloseq package or the MicrobiomeAnalyst tool (https://www.microbiomeanalyst.ca) (25).

Differentially abundant ASVs between non-SS and SS subjects were identified using the following two methods: DESeq2 (26) and metagenomeSeq (27). DESeq2 (v.1.33) and metagenomeSeq (v.1.35) were performed in the R environment. Results from different methods were collected, and intersections at the genus level were plotted in a Venn diagram using InteractiVenn (28) (http://www.interactivenn.net). Logistic regression was performed using the run_sl function in the MicrobiomeMarker package (29). Briefly, the abundance of each genus was transformed to log10(1+x), and the relative abundance was normalized with total sum scaling. The top eight genera, which had a feature importance score of more than 10, were used to perform logistic regression. Logistic regression performed with more than eight genera did not significantly improve analysis outcomes. The receiving operating characteristic curve was plotted with the plot_sl_roc function in the MicrobiomeMarker package.

The 16S rRNA sequencing data were deposited in the Sequence Read Archive of the National Center for Biotechnology Information with the BioProject accession number of PRJNA804331.

All study patients were female with a mean age of 51.4 years. The demographic and clinical characteristics of the study patients are shown in Tables 1 and 2 . The SS sicca and non-SS sicca groups comprised 23 and nine patients, respectively ( Table 1 ). Among the SS subjects, 15 (65.2%) tested positive for anti-SSA/Ro antibodies and 23 (100%) exhibited lymphocyte infiltration in the LSGB samples. Of the 23 SS subjects, 18 (78.3%) and 13 (56.5%) exhibited dry mouth and eyes, respectively. All patients in the non-SS sicca group exhibited dry mouth, while 5 (55.6%) exhibited dry eyes. All patients in the SS cohort group were confirmed to show more than a score of 4 according to the weight of the new criteria established by the ACR and EULAR in 2016, whereas all non-SS cohort patients showed negative findings in both anti-SSA/Ro antibodies and lymphocyte infiltration in the LSGB samples, with a total score less than 4.

UWMSFR in the SS (0.28 ± 0.29 mL/min) group was not significantly different from that in the non-SS group (1.03 ± 0.79 mL/min) ( Table 2 ). Salivary gland scintigraphy parameters, such as UR, MA, and Tmin were not significantly different between the two groups. However, MS was significantly reduced in the SS cohort. Next, the MR sialography features of PGs were analyzed ( Figures 1A – E and Table 2 ) and radiological structural deformities were compared between the two groups. In the PGs of SS patients, the sialectactic dilated duct was prominent in 12 (63.16%) glands. However, none of the non-SS subjects exhibited PG sialectasis. The fat deposition in the submandibular gland (SMG) was increased in the SS cohort, whereas PG fat depositions were not significantly different between the two groups ( Figures 1F – I and Table 2 ).

To investigate the correlation between the microbiome and the biological features of SS, PG saliva collected through PG lavage was subjected to bacterial 16S rRNA sequencing. The observed ASV, which was inferred using DADA2, in the SS group was higher than that in the non-SS group. Alpha diversity indices, including Shannon and Simpson indices, were not significantly different between the two groups ( Figure 2A ). Meanwhile, alpha diversity was significantly downregulated in the anti-SSA-positive group. This indicated that the microbiota in PG saliva exhibited low diversity and that it was correlated with anti-SSA-positivity ( Figure 2B ). The observed ASV and alpha diversity were not significantly different when the groups were compared based on the focus score ( Supplementary Figure 1A ).

Next, the microbial abundance at the phylum level was examined in patients with SS. The abundance of the phylum Firmicutes in the anti-SSA-positive group was higher than that in the anti-SSA-negative group. In contrast, the abundance of the phyla Actinobacteria and Proteobacteria in the anti-SSA-positive group was significantly lower than that in the anti-SSA-negative group ( Figure 2C ). In addition, the abundance of the phylum Firmicutes in patients with SS was also higher than that in non-SS subjects ( Supplementary Figure 1B ). The observed ASV was significantly correlated with the abundance of the phylum Bacteroidetes ( Figure 2D ). Meanwhile, the abundance of the phylum Firmicutes was negatively correlated with the Shannon index, whereas that of the phyla Bacteroidetes and Actinobacteria was positively correlated ( Figure 2E and Supplementary Figure 1C ). Additionally, ASV and microbial diversity were negatively correlated with log-scaled Firmicutes/Bacteroidetes ratio (F/B ratio, Figures 2D, E ). Rarefaction curves further indicated that our sequencing depth was sufficient to reach maximum coverage ( Supplementary Figure 1D ). Decreased diversity and alteration of phylum abundance in anti-SSA-positive group were also observed when only SS samples were used for analysis ( Supplementary Figures 2A – D ). Altogether, these findings suggest that the PG saliva microbiome exhibits low diversity dysbiosis in patients with SS and that this is correlated with biological markers of SS (anti-SSA antibodies).

To evaluate the clinical implications of the PG saliva microbiome in patients with SS, the correlation of the microbiota composition with clinical (symptoms), functional (SFR), and radiological (MR sialography) features was examined. The ASV and alpha diversity were not significantly different between the groups based on clinical symptoms, including dry eye, dry mouth, and parotid swelling or pain ( Figures 3A – C ). Among patients who were tested, SWMSFR was positively correlated with the abundance of the phylum Bacteroidetes ( Figure 3D ). Interestingly, we found that a significant decrease in the Shannon index (less than 3) was only observed in patients with sialectactic duct ( Figure 3E ). The Shannon index, the abundance of Firmicutes, and F/B ratio were significantly different when the groups were compared based on PG ductal types 0 (healthy) – 1 (stenosis) versus 2 (sialectasis) ( Figure 3F ). Additionally, fat deposition in the left or right PG was positively correlated with the Shannon index ( Figure 3G ). Other functional and radiological features, such as salivary gland scintigraphy parameters (US, MA, MS, and Tmin), PG sialectasis grade, PG fat stage, and SMG fat stage were not correlated with the microbiota composition in patients with SS ( Supplementary Figures 3A, B , and data not shown). These findings indicate that distinct clinical features, including PG ductal deformity (sialectasis) and consequent salivary stasis, may be correlated with PG microbial dysbiosis and enrichment of Firmicutes in patients with SS. Additionally, high SFR without ductal dilation and fatty degeneration may be correlated with the relative abundance of Bacteroidetes and microbial diversity, respectively.

Next, the ability of the abundance of microbial genera to distinguish between non-SS and SS subjects was examined. The principal coordinates analysis plot indicated that disease status contributed to 5.9% of PG saliva microbiota composition variation with statistical significance ( Figure 4A ). The microbiome in the non-SS group mainly comprised relatively diverse genera, such as Brevundimonas, Streptococcus, Haemophilus, and Lactobacillus ( Figure 4B ), whereas that in the SS group was characterized by increased abundance of Streptococcus and decreased microbial diversity ( Figure 4C ). To identify potential biomarkers to distinguish between SS and non-SS subjects, the following two different differential analysis methods were used: DESeq2 and metagenomeSeq. RNA sequencing-based DESeq2 analysis revealed that several genera exhibited differentially abundance between the non-SS and SS groups ( Figure 4D ). Furthermore, the metagenomeSeq-driven volcano plot revealed that the abundance of Bacteroides vulgatus, Staphylococcus aureus, Methylorubrum extorquens, Blautia coccoides, and Lactobacillus johnsonnii was significantly downregulated in the SS group ( Figure 4E ). Although a biomarker that was specifically enriched in the SS group was not identified, Lactobacillus was the common genus identified in DESeq2 and metagenomeSeq analyses. This suggested that the abundance of Lactobacillus is downregulated in SS subjects but not in non-SS subjects ( Figure 4F ).

Furthermore, the logistic regression model was used to determine whether the abundance of the microbiome at the genus level could be used to predict SS. The model was constructed with the top eight genera with the area under the curve value of 0.79 ( Figure 5A ). Among the tested genera, four genera, including Streptococcus, exhibited featured importance scores in patients with SS. Similarly, four genera, including Lactobacillus, exhibited featured importance scores in non-SS subjects ( Figure 5B ). The relative abundance of each genus also exhibited a similar enrichment pattern based on the feature importance score ( Figure 5C ). This suggested that these genera can differentiate SS subjects from non-SS subjects.

A large overlap in the microbial composition was observed between the non-SS and SS groups, while separate microbial clusters were observed in some SS subjects ( Figure 6A ). Thus, this cluster was termed the “dysbiosis cluster” and used for further analysis. The microbiome in the dysbiosis cluster was highly enriched with the genus Streptococcus, while various taxa were enriched in the non-dysbiosis cluster ( Figures 6B, C ). Differential analysis of DESeq2 and metagenomeSeq results revealed that the abundance of Streptococcus oralis was high, while the abundance of other species, including Bacteroides vulgatus, Bifidobacterium longum, and Faecalibacterium prausnitzii, was markedly downregulated in the dysbiosis cluster ( Figures 6D, E ). These findings indicate that Streptococcus and Lactobacillus are potential biomarkers for SS and non-SS subjects, respectively, at the genus level.