Patients and controls were age and gender-matched and randomly allocated across the different experiments. All pSS patients fulfilled the AECG classification-criteria for pSS (30). Patients that did not fulfill the pSS classification-criteria, but presented with dryness-complaints without a known cause, in the absence of any generalized autoimmune disease were classified as non-Sjögren’s sicca (nSS) patients ( Table 1 ). The study was approved by the medical ethics committee of the University Medical Center Utrecht (METC no. 13-697). All patients gave their written informed consent in accordance with the declaration of Helsinki.

Monocytes were isolated from peripheral blood mononuclear cells by magnetic-activated cell sorting using the CD14⁺ isolation kit (Miltenyi Biotec) according to the manufacturer’s instructions. The frequency of monocytes subsets and the purity of the isolated monocytes was measured by flow-cytometry to insure consistent results ( Supplementary Table 1 ). Monocytes subsets were defined as the following: classical monocytes (CD14⁺CD16^-), intermediate monocytes (CD14⁺CD16⁺), and non-classical monocytes (CD14^-CD16⁺) ( Supplementary Figure 1 ). The purity of the isolated samples was (median [range]) 96% [89–99%], and there were no significant differences in cell purity between any of the groups. Cells were lysed in RLTplus buffer (Qiagen) supplemented with 1% beta-mercaptoethanol and total RNA was purified using AllPrep Universal Kit (Qiagen), according to the manufacturer’s instructions. RNA concentration was assessed with Qubit RNA Kit (Thermo Fisher Scientific) and RNA integrity was measured by capillary electrophoresis using the RNA 6000 Nano Kit (Agint Technologies); all samples had RIN-score >7.0.

RNA sequencing was performed at the Beijing genomics institute, using an Illumina HiSeq 4000 sequencer (Illumina) following the standard manufacturer’s protocols. About 20 million 100bp paired-end reads were generated for each sample. RNA‐sequencing data analysis was performed as previously described (31).

Briefly, quality control of the reads was assessed with FastQC tool (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and all samples passed the quality check. Next, reads were aligned to the human genome using STAR aligner (32). The Python package HTSeq (33) was used to calculate the read counts for each annotated gene. Differentially expressed genes (DEGs) were identified using Bioconductor/R package DESeq2 (34); Wald’s test was used to identify DEGs in each pair-wise comparison performed between the three groups (HC, nSS, and pSS) and likelihood ratio test to identify DEGs considering multiple groups. P-values were corrected using a false discovery rate (FDR) of 5% according to the Benjamini and Hochberg method. Differences in gene expression with a corrected p-value<0.05 were considered differentially expressed. Variance stabilizing transformation was applied to obtain normalized gene counts (variance stabilized data), which were used for subsequent analyses.

Gene set enrichment analysis was conducted with Broad Institute software (35), to explore the contribution of each monocyte subset to the transcriptomic profile identified in monocytes from pSS patients. Thousand random permutations were performed to calculate the normalized enrichment score and FDR-corrected p-value. Gene sets were considered significantly enriched with a corrected p-value<0.05. Gene sets were obtained from a publicly available microarray dataset which identified the representative gene signature of each monocyte subsets (classical, intermediate, and non-classical monocytes) (10).

To study the inter-dependence between genes and their interactions, a gene co-expression network was performed using weighted gene co-expression network analysis (WGCNA) Bioconductor/R package (36). Gene co-expression networks were constructed using all differentially expressed genes (corrected p<0.05) in both comparisons, nSS or pSS vs. HC ( Supplementary Figure 2 ). Networks were defined by first performing an unsigned pairwise spearman’s correlation between genes and subsequently transforming the co-expression similarities into an adjacency matrix and scaling the adjacency matrix to achieve a scale-free topology (scaling power = 9, selected to have a network that fits the scale-free topology criterion). The identified networks were labelled using a conventional color scheme (blue, brown, turquoise and yellow) and further designated as gene signatures.

To understand the biological meaning of the identified signatures, pathway enrichment analysis was performed in the different sets of genes using the compareCluster function in the ReactomePA package (37) to compare and plot the pathways enriched in the different signatures.

For each signature, the expression profile was summarized by the first principal component of the signature expression level to calculate the module eigengene expression, which represents the average expression level of the genes within the signature. Next, to identify signatures that were correlated with relevant clinical features of pSS patients, a module-trait association was calculated using the pearson’s correlation coefficient between the module eigengene and the clinical traits.

To gain further insight into the signatures and the genes that have the closest relationship with the clinical features, an intramodular analysis was performed. Module membership defined as the degree of correlation between the gene and the signature, and gene significance calculated as the absolute value of the association between the gene expression profile and each clinical trait were used. Correlation with a r value higher than 0.3 and a p-value<0.05 were plotted and selected for hub-gene identification. Hub-genes were identified in the brown, yellow and turquoise signature and defined by a module membership higher than 0.8 and a gene significance higher than 0.4. Next, hierarchical cluster analysis of the hub-gene expression based on Ward’s method with Euclidian distance was performed using the MetaboAnalyst5.0 online software (http://www.metaboanalyst.ca/) to identify patients with a similar inflammatory profile.

The IFN-score was calculated using the gene expression of five IFN-induced genes (IFI44L, IFI44, IFIT3, LY6E, and MX1), after z-score normalization, as previously described (22).

As functional annotation of the turquoise signature revealed enrichment in genes involved in TLR-signaling we set out to confirm these alterations on a functional level. For this purpose, whole blood was diluted 1:1 in RPMI-1640 medium with 1% L-glutamine (both from Thermo Fisher Scientific) and stimulated with different TLR ligands ( Supplementary Table 2 ). 1h after stimulation, 10μg/mL of Brefeldin A (Sigma) was added. After a total of 6h of stimulation cells were stained with the extracellular antibodies and after washing, fixation and permeabilization with FIX&PERM (Thermo Fisher Scientific), cells were stained with anti-TNF-α ( Supplementary Table 1 ). Data acquisition was performed using a BD LSRFortessa (BD Biosciences) and data were analyzed using FlowJo software (Tree Star).

Monocytes from HC were isolated as previously described, and treated with serum of HC or pSS patients (20% v/v) for 20h, either alone or in the presence of TLR4 ligand (0.1μg/mL of LPS; Invivogen) for the last 17h of culture. Simultaneously, monocytes were pre-incubated for 1h at 37°C with 5μg/ml of either anti-IFNα/βR2 (clone MMHAR-2; PBL Assay Science) or its respective isotype control (IgG2a, Thermo Fisher Scientific), and then treated with serum of pSS patients (20% v/v) for 24h. Supernatants were harvested and cells lysed in RLTplus buffer (Qiagen). RNA isolation was performed as previously described. Quantitative-PCR reactions were performed using SYBR Select Master Mix (Applied Biosystems) on the Quantstudio 12k system (Thermo Fisher Scientific). Relative gene expression levels were normalized to the geometric mean expression of two housekeeping genes: B2M and ACTB ( Supplementary Table 3 ). The relative fold change (FC) of each sample was calculated in relation to the average ΔCt of the HC group or to the ΔCt of the respective isotype condition (reference), according to the formula FC = 2^−ΔΔCt, where ΔΔCt = ΔCt sample —ΔCt reference. TNF-α levels were measured in the supernatants by ELISA (Diaclone), following the manufacturer’s instructions.

Differences in gene expression with a corrected p-value<0.05 were considered statistically significant. Kruskal-Wallis H test was used to analyze the differences between the nSS, pSS, and HC groups and between patient clusters. Fisher’s exact test was used to compare categorical variables. Differences between pSS patients and HC were assessed using Mann-Whitney U-test. Changes observed after IFNα/βR2 blocking were analyzed using Wilcoxon matched-pairs signed rank test. Statistical analyses and data visualization were performed using Python and R, Graphpad Prism (GraphPad Software), and MetaboAnalyst 5.0 (38). Differences were considered to be statistically significant at p<0.05.

To better understand the contribution of monocytes to the development of pSS and to identify monocyte-specific processes dysregulated in pSS, we performed RNA-sequencing profile of circulating monocytes of pSS and nSS patients compared to HC ( Table 1 ). Differential expression analysis revealed that among the studied groups, 6.036 genes were differentially expressed with a nominal p-value ≤ 0.05, of which 3.126 were upregulated and 2.910 downregulated. The majority of the differentially expressed genes (3.408 DEGs) were found in the pSS vs. HC comparison, representing 56% of all the DEGs, and out of these 1.731 were uniquely differentially between pSS and HC. In addition, 29% off all DEGs (1.757 genes) were differentially expressed in nSS patients vs. HC, and only 465 (26%) were exclusively found differentially expressed between nSS and HC. Furthermore, in the comparison between pSS and nSS, out of the 871 DEGs only a small number of genes (280/32%) was uniquely differential between pSS vs. nSS ( Figure 1A ). Overall, the transcriptomic profile of pSS monocytes revealed clearly distinct expression profiles compared to HC monocytes, but relatively similar profiles to nSS monocytes. To identify the most robust and consistently altered genes, those differentially expressed with an FDR corrected p-value ≤ 0.05 and with at least moderately high expression levels (base mean expression >100) were selected for further analysis ( Supplementary Figure 2 ).

To gain further insights into the molecular signature of pSS monocytes, we performed gene set enrichment analysis to evaluate the contribution of each monocyte subset to the transcriptional profile found in pSS monocytes. To test this, we compared the DEGs in pSS patients that met the selection criteria ( Figure 1B ) with publicly available gene expression profiles of the three different monocyte subsets (10). Enrichment analysis revealed that genes highly expressed in intermediate and non-classical monocytes were significantly enriched in the transcriptional profile of pSS monocytes ( Figure 1C ). In line with this, flow-cytometry analysis of the monocyte subsets confirmed that the non-classical monocyte subset was increased in the circulation of pSS patients when compared with nSS patients and HC ( Figure 1D ). Thus, the transcriptomic alterations found in pSS monocytes are consistent with an increased frequency of non-classical monocytes.

To further analyze the functional consequences of the transcriptomic changes in pSS monocytes, we used weighted gene co-expression network analysis (WGCNA) to construct gene-correlation networks and identify signatures of co-expressed genes with pSS-associated biological functions. Network analysis using the genes that were differential between any of the three groups compared ( Supplementary Figure 2 ) identified 4 gene signatures: blue (358 genes), brown (239 genes), turquoise (737 genes), and yellow (148 genes), each with distinct expression patterns ( Figure 2A and Supplementary Table 4 ). 18 of the DEGs were not allocated to any of the previous signatures and therefore excluded from further analysis.

Functional annotation indicated that the blue signature was associated with translation mechanisms, the brown signature with IFN-signaling, and the turquoise signature with TLR signaling ( Figure 2B ). No pathways were enriched for the yellow signature. Next, we further characterized the identified signatures, by assessing the module eigengene expression, given by the first principal component of the signature expression level, in HC, nSS, and pSS patients. Overall, the expression of genes from the blue signature was decreased in nSS and pSS patients and contains a large number of ribosomal protein genes, and a range of initiation factors and RNA polymerase subunit genes. Conversely, the expression of the brown, turquoise, and yellow signatures was largely increased in both patient groups, represented by an increased eigengene expression. The expression patterns of all four gene signatures in the nSS group were in general intermediate between the HC and pSS groups ( Figure 2C ).

To better understand the biological relevance of the different signatures for the disease, we next assessed the correlation between clinical traits ( Table 1 ) and the eigengene expression in each signature. Eigengene expression is an ideal parameter to correlate with the clinical traits, since it represents the first principal component of the co-expressed gene network and thus accounts for most of the variance in gene expression within the signature. Genes from the brown signature had the most correlations with the clinical traits, including the strongest positive correlation with ESR (r = 0.51, p = 0.003), and with ESSDAI (r = 0.37, p = 0.042). In addition, the lymphocytic focus score (LFS), which reflects salivary gland inflammation, moderately correlated with the turquoise and yellow signatures (r = 0.43, p = 0.015 and r = 0.40, p = 0.024, respectively), whereas only the yellow signature moderately correlated with Schirmer test, a measure for lacrimal gland function (r = 0.44, p = 0.013). Of note, moderate negative correlations were found with C4 levels for the brown and yellow signatures (r = -0.51, p = 0.003 and r = -0.43, p = 0.016, respectively) ( Figure 2D ).

Next, we performed an intramodular analysis to seek out the biological significance of each gene within the respective signature for the clinical features that correlated with the signature expression ( Figure 2D ). The genes within the brown, and turquoise signatures significantly correlated with sIgG, ESR, and complement C3 levels, parameters usually associated with systemic inflammation. Likewise, the genes within the yellow signature correlated with complement C4 levels and with Schirmer test, indicative of lacrimal gland dysfunction ( Figure 2E ). No correlations were observed for LFS and ESSDAI. Genes with a module membership (a higher number indicates stronger correlation of the gene with the other genes in the signature) higher than 0.8 and a gene significance (a higher number indicates stronger correlation with the selected clinical trait) higher than 0.4 were selected as hub-genes for each signature independently. Analysis of the brown signature revealed 43 hub-genes with gene significance for sIgG, and 68 hub-genes with gene significance for ESR. Out of these, 31 hub-genes were found with gene significance for both parameters ( Figure 2E and Supplementary Table 5 ). Within the yellow signature, 16 hub-genes were identified with gene significance for Schirmer test and 13 were found with gene significance for C4 levels, of these 7 hub-genes were identified with gene significance for both the Schirmer test and C4 levels ( Figure 2E and Supplementary Table 6 ). In the turquoise signature, analysis of gene significance for C3 complement levels yielded 107 hub-genes ( Figure 2E and Supplementary Table 7 ). Thus, the molecular signatures of circulating monocytes from pSS patients seems to reflect the ongoing systemic inflammation as result of local damage.

Given the fact that hub-genes are highly associated with the clinical features of pSS, we next investigated whether the hub-gene expression could be used to identify patients with a similar inflammatory profile. Unsupervised hierarchical cluster of the selected hub-genes ( Supplementary Tables 5 – 7 ) identified four distinct clusters of nSS and pSS patients. Cluster 1 mainly included patients diagnosed as pSS patients (5 out of 6), and all patients were characterized by an increased expression of the brown signature ( Figure 2F and Supplementary Figure 3A ). As a matter of fact, patients in cluster 1 showed the highest IFN-score compared to the patients included in other clusters ( Figure 2G ). Clustering analysis identified a group of patients included in cluster 2 with a pronounced downregulation of the yellow signature hub-genes, and a trend towards a lower LFS, ESR, sIgG and IFN score ( Figure 2G and Supplementary Figure 3A ). Patients in cluster 3 were characterized by an overexpression of the turquoise signature and displayed an intermediate profile for ESR, IFN-score and frequency of anti-SSA antibodies, as compared with patients in cluster 1 and 2 ( Figure 2G and Supplementary Figure 3A ). Patients in cluster 4 did not display a specific signature profile ( Supplementary Figures 3A, B ). Furthermore, comparison of systemic autoimmune and inflammatory parameters between clusters showed that patients in cluster 1 have an increased frequency of anti-SSA antibodies and presented with higher ESR. Moreover, a trend towards a higher sIgG and lower C3 and C4 levels was observed in these patients compared to the patients in the other clusters ( Figure 2G and Supplementary Figure 3B ). Finally, patients in cluster 1 have higher frequency of non-classical monocytes compared to the patients in cluster 3 and 4 ( Supplementary Figure 3C ).

Considering that the functional annotation of the turquoise signature revealed enrichment in genes involved in TLR-signaling, we next assessed whether the altered transcriptional profile of pSS-monocytes could impact monocyte activation. To confirm this, whole-blood from pSS patients and HC was stimulated with a broad panel of toll-like receptor ligands and the intracellular expression of TNF-α was evaluated by flow cytometry ( Supplementary Figure 4 ). Interestingly, the intracellular TNF-α expression in monocytes of pSS patients, as determined by the median fluorescence intensity (MFI), was significantly increased upon activation with TLR4, TLR5, and TLR7/8 ligands compared to HC monocytes ( Figure 3A ). We did not observe differences in TNF-α expression between HC and pSS monocytes upon TLR2/1, TLR2/6, and TRL3 activation ( Figure 3A ).

As the molecular signatures of monocytes from pSS patients and the expression of their hub-genes were mainly associated with systemic inflammatory markers, we sought to confirm whether the transcriptomic changes found in pSS-monocytes could be induced in monocytes from healthy individuals by sera of pSS patients. To this end, we selected hallmark genes to represent each signature, and ultimately the transcriptomic profile of pSS monocytes, by selecting those hub-genes with the highest degree of correlation with each signature (high module membership), the highest biological relevance (high gene significance) and highest fold-change of expression differences between pSS and HC ( Supplementary Table 8 ). Treatment with serum from pSS patients upregulated both hallmark genes of the brown signature (MX1, IFITM1), two of the three turquoise hallmark genes (TNFSF10, IRF2), one of the two yellow signature hallmark genes (IRF9), and downregulated one of the three blue signature hallmark genes (RPL5) ( Figure 3B ). In addition, we tested whether the inflammatory mediators present in the serum of pSS patients could prime healthy monocytes for higher TNF-α production upon TLR activation. Due to limited amount of serum, we were only able to test this upon TLR4 activation. We observed that the changes induced by serum from pSS patients led to a significantly higher TNF-α production upon TLR4 stimulation ( Figure 3C ).

As type-I IFNs are considered to be crucial mediators of pSS pathogenesis (39), we tested whether the transcriptomic alterations were mediated by type I IFNs. Interferon-α/β receptor (IFNAR) blockade completely abrogated the pSS-serum effect on the brown gene signature, and also reverted the pSS-serum effects on TNFSF10 and IRF9 expression, hallmark genes of the turquoise and yellow signature, respectively. In addition, RPL5 expression, hallmark of the blue signature, was upregulated upon IFNAR blockade, although not statistically significant. Even though the expression of RPL15 and EEF1B2, hallmarks of the blue signature, were not significantly downregulated upon treatment with serum from pSS patients, IFNAR blockade significantly increased the expression of both genes. ( Figure 3D ). Thus, type-I IFNs importantly drive the observed effects of pSS serum on the identified signature hallmark genes in monocytes.