The Transcriptome of Paired Major and Minor Salivary Gland Tissue in Patients With Primary Sjögren’s Syndrome

Background While all salivary glands (SGs) can be involved in primary Sjögren’s syndrome (pSS), their respective role in pathogenesis remains unclear. Our objective was to assess immunopathway activation in paired parotid and labial gland tissue from biopsy-positive and biopsy-negative pSS and non-SS sicca patients. Methods Paraffin-embedded, paired parotid and labial salivary gland tissue and peripheral blood mononuclear cells were obtained from 39 pSS and 20 non-SS sicca patients. RNA was extracted, complementary DNA libraries were prepared and sequenced. For analysis of differentially expressed genes (DEGs), patients were subdivided based on fulfillment of ACR-EULAR criteria and histopathology. Results With principal component analysis, only biopsy-positive pSS could be separated from non-SS sicca patients based on SG gene expression. When comparing the transcriptome of biopsy-positive pSS and biopsy-negative non-SS sicca patients, 1235 and 624 DEGs (FDR<0.05, log2FC<-1 or >1) were identified for parotid and labial glands, respectively. The number of DEGs between biopsy-negative pSS and non-SS sicca patients was scarce. Overall, transcript expression levels correlated strongly between parotid and labial glands (R2 = 0.86, p-value<0.0001). Gene signatures present in both glands of biopsy-positive pSS patients included IFN-α signaling, IL-12/IL-18 signaling, CD3/CD28 T-cell activation, CD40 signaling in B-cells, DN2 B-cells, and FcRL4+ B-cells. Signature scores varied considerably amongst pSS patients. Conclusion Transcriptomes of paired major and minor SGs in pSS were overall comparable, although significant inter-individual heterogeneity in immunopathway activation existed. The SG transcriptome of biopsy-negative pSS was indistinguishable from non-SS sicca patients. Different patterns of SG immunopathway activation in pSS argue for personalized treatment approaches.


INTRODUCTION
In primary Sjögren's syndrome (pSS), salivary and lacrimal glands are considered main targets of the disease (1). A pathogenic hallmark of pSS is B-cell hyperactivity, reflected by the presence of autoantibodies (anti-SSA/Ro, anti-SSB/La, rheumatoid factor), hypergammaglobulinemia, and the increased risk of developing B-cell lymphoma (2). In addition to the presence of disease-associated autoantibodies, another important diagnostic feature of pSS is periductal focal lymphocytic sialadenitis (FLS) in salivary glands (SGs) (3). FLS with a focus score (FS) ≥1 indicates a positive SG biopsy (4). Other histologic features associated with pSS are the presence of lymphoepithelial lesions (LELs) and a plasma cell shift (relative decrease of IgA plasma cells with influx of IgM and/or IgG plasma cells). Furthermore, within the SG infiltrates, ectopic lymphoid structures can develop that may contain germinal centers (5). While the intimate interaction between glandular epithelium and immune cells is likely critical to pSS pathogenesis (6), a clear correlation between immune cell infiltration and main patient-reported symptoms (i.e. dryness, pain and fatigue) or salivary flow rates is lacking (7)(8)(9). There is, however, a positive association between FS and the amount of atrophy and fibrosis (8,10), indicating that infiltrates can lead to glandular tissue damage.
SG infiltrates are mainly composed of CD4 + T-cells and Bcells, but may also harbor other immune cells, such as dendritic cells and natural killer cells (11). T-cell dependent B-cell hyperactivity and type-I interferon (IFN) activation are both thought to play central roles in the immunopathogenesis of pSS (2,12,13). Peripheral B-cell and type-I IFN activation signatures, measured by whole blood transcriptomics, have been associated with separate clinical phenotypes (14). Other ways to stratify pSS patients based on immune signatures in blood have also been proposed (15,16). Alternatively, patient stratification could be based on immune signatures in SG tissue, although so far only smaller-sized transcriptomic studies with labial salivary gland (LSG) tissue have been performed (17)(18)(19). While all types of SGs can be involved in pSS, LSGs are most commonly biopsied for diagnostic purposes, with parotid salivary glands (PSGs) as a suitable alternative (20). The respective roles and contributions of different glands to immunopathogenesis remain unclear. Such knowledge would advance our understanding of the disease process and is essential for patient stratification and clinical trial design with targeted immunotherapies. We therefore analyzed the transcriptome of paired PSG and LSG tissue, as well as matched peripheral blood mononuclear cells (PBMCs) from pSS and non-SS sicca patients. We compared immune signatures between biopsy-positive and -negative patient groups and tissues and correlated our findings to clinical and histopathological parameters.

Study Population
We included 39 consecutive patients with pSS who participated in our inception cohort and fulfilled 2016 ACR-EULAR criteria (21). Twenty non-SS sicca patients from the same cohort, age-and sexmatched to the pSS group, were also included. Patients agreed to undergo both a parotid and labial gland biopsy. Matched cryopreserved PBMC, whole blood, and serum samples, collected at inclusion, were also available. Paired parotid and labial gland samples were always collected on the same day. Blood samples were collected when patients visited the rheumatologist around the same time (maximum of one month before or three months after the biopsy). No clinical interventions were initiated before biopsies were taken. Informed consent was obtained according to the Declaration of Helsinki. The study was approved by the Medical Research Ethics Committee of the University Medical Center Groningen (METc2013.066).

Subgroup Analysis Based on Histopathology
We divided patients into four groups based on fulfillment of 2016 ACR-EULAR criteria and histopathological phenotype. The rational for this grouping was based on observations that non-SS sicca patients may have mild lymphocytic infiltration and pSS patients may have a negative biopsy result (i.e., focus score<1.0). For the purpose of this study, a positive biopsy was defined as FLS with focus score≥1 and/or presence of LELs, because LELs are highly specific for pSS compared with non-SS sicca (22). Biopsies with sclerosing chronic sialadenitis, mucosa-associated lymphoid tissue (MALT) lymphoma, or severe atrophy and PSG biopsies containing a lymph node were excluded.

RNA Extraction and Sequencing Library Preparation
RNA was extracted from 11 formalin-fixed (4%), paraffinembedded tissue sections (4mm) per biopsy. After scraping the tissue from the slides and deparaffinization, tissue was digested in buffer PKD (Qiagen, USA). RNA was extracted using the RNAeasy FFPE kit, according to manufacturer's protocol (Qiagen, USA). For RNA extraction from PBMCs, 3x10^6 cells were thawed and rested for 1h in RPMI medium supplemented with 10%FBS and Antibiotic-Antimycotic (all from Gibco, UK) at 37°C. After resting, cell lysis and RNA isolation were performed using the RNAeasy Mini kit, according to manufacturer's protocol (Qiagen, USA). RNA was immediately stored at -80°C until further use. Complementary DNA library preparation was performed by using TruSeq Stranded Total RNA Library Prep Gold (Illumina, USA, Cat#20020598), following manufacturer's recommendations. Prepared libraries were sequenced on a HiSeq 2500 System (Illumina, USA).

Transcriptomic Data Analysis
Data quality assessment and principal component analysis (PCA) were performed to understand the main source of variability. Differential expression analysis was performed using the limma package (v3.42.2) in R (v3.6.1) to identify differentially expressed genes (DEGs) across subgroups. Adjusted P-values were calculated using the Benjamini-Hochberg method; a gene was considered differentially expressed if the absolute fold change was >2 or <0.5 and the false discovery rate-adjusted P-value was <0.05. The MetaCore from Clarivate Analytics pathway database was used for pathway enrichment analysis. Additionally, gene signatures were established based on either published data (i.e., DN2 B-cell signature, FcRL4+ B-cell signature (23,24)) or in vitro modulation of healthy human PBMCs or purified cells. Details on gene signature generation and composition are provided in Supplementary Tables 1 and 2. The single-sample Gene Set Enrichment Analysis algorithm was used to calculate the score of pre-defined gene signatures per sample (25). For evaluation of the correlation between average gene expression in the two types of salivary gland, only patients with paired samples were included.

Peripheral Blood Cytokine Analysis
Concentrations of 84 cytokines in serum were assessed using multiplexed bead-based immunoassays (Merck Millipore, USA). Concentrations of myxovirus resistance protein A (MxA), a surrogate marker for type-I IFN activity (26), were measured in lysed whole blood using an in-house ELISA (see Supplementary Methods).

Histological Assessments
For each biopsy, the focus score and area fraction of the infiltrate were calculated. To define the SG area occupied by the infiltrate, sections were immunohistochemically stained with anti-CD45 (Dako, clone 2B11+PD7/26). The relative amount of CD45positive infiltrate was assessed in relation to the total amount of tissue (excluding larger areas of fat and fibrotic tissue) by a pixel-based digital image analysis in QuPath (27). Table 1. Out of the 59 patients included in this study, one patient had withdrawn consent for the labial gland biopsy after inclusion.

Patient characteristics are depicted in
Of the remaining 58 patients with a paired PSG and LSG biopsy, we excluded seven SG samples (each from individual patients) for the following reasons: confirmed MALT lymphoma (n=1; parotid), sclerosing chronic sialadenitis (n=1; labial), presence of a lymph node (n=3; parotid), or severe atrophy (n=2; parotid) from the analysis. Thus, in total 52 PSG and 57 LSG samples were included in the analyses with 51 paired biopsies of 34 pSS and 17 non-SS sicca patients. Numbers of included samples per subgroup are shown in Table 2. PBMC, serum, and whole blood samples were available from 59, 57, and 52 patients, respectively. Because the biopsy-positive non-SS group, as expected, contained only few samples, this group was excluded from statistical analyses.

Only Biopsy-Positive pSS Patients Show an Altered SG Transcriptome
First, we performed a PCA to explore whether PSG and LSG samples from pSS and non-SS sicca patients could be distinguished by their transcriptome. PCA revealed that only pSS patients with a positive biopsy (group-IV) could be separated from biopsy-negative pSS (group-III) and non-SS sicca patients (groups I&II) based on gene expression analysis in either PSG or LSG tissue ( Figure 1A). Distinct clustering was associated with a relative increase in immune infiltrate size, measured by the CD45 + area of parenchymal tissue ( Figure 1A).
Based on PCA results, we then compared the transcriptome of biopsy-positive pSS (group-IV) with biopsy-negative non-SS (group-I) to look for genes that were differentially expressed. For the PSG, 1041 up-regulated (FDR<0.05, log2FC>1) and 194 down-regulated (FDR<0.05, log2FC<-1) genes were identified. For the LSG, these numbers were 581 and 43, respectively. The top 20 up-regulated genes in both tissues were mostly B-cell or Tcell related (full list of DEGs in Supplementary Tables 3 and 4). Generally, large overlap between upregulated DEGs in PSGs and LSGs was observed ( Figures 1B, C). There was a strong correlation between average gene expression in PSGs and LSGs (R 2 = 0.86, Figure 1C) and between transcript expression fold changes (R 2 = 0.59) when taking all genes into account. Notably, also the focus score and area fraction of CD45+ cells showed a good correlation between paired PSG and LSG samples from pSS patients (Spearman's rho=0.65 and 0.66, respectively; p<0.001; Supplementary Figure 1), although higher values were observed in LSG vs. PSG tissue ( Table 1). Despite a higher degree of infiltration in LSG tissue, we observed a tendency towards higher transcript fold changes in PSG tissue when comparing biopsypositive pSS and biopsy-negative non-SS patients (Supplementary Figure 2). Genes with poorly correlated expression levels between the two SG types were not associated with immune function and mostly related to SG biology (e.g., genes encoding cystatins, ATPases, carbonic anhydrase VI, and a-amylase; data not shown).
Secondly, we compared the transcriptome of biopsy-positive pSS (group-IV) and biopsy-negative pSS (group-III) patients. DEGs were very similar to those identified in the previous comparison (group-IV vs. group-I; data not shown). In marked contrast, when we compared biopsy-negative pSS (group-III) with biopsy-negative non-SS sicca (group-I) patients, only two genes in LSG and none in PSG tissue were differentially regulated ( Figure 1D). While transcriptomes of these groups were almost indistinguishable, at the functional level we noted some difference, since unstimulated whole salivary flow (UWSF) rates tended to be higher in group-I vs. group-III for the LSG (not PSG) comparison (p=0.13). There were no differences in stimulated WSF (SWSF) rates between the groups (Supplementary Table 5).

A Key Set of Immunopathways Is Active in Both Glands
Next, we analyzed whether specific pathways were enriched amongst up-regulated genes in SGs of biopsy-positive pSS compared to non-SS sicca patients. A pathway database was used for pathway enrichment analysis of DEGs. As expected, the majority of upregulated pathways involved T-cells or B-cells ( Figure 2A). Identified pathways were largely similar between PSG and LSG analyses. A schematic overview of upregulated genes with known immune function is provided as Figure 2B.
Additionally, to obtain individual sample data on relevant biological processes, gene signature analysis was also performed. These signatures were selected based on current knowledge of immune pathways involved in pSS pathogenesis (12,13,28,29). For both PSG and LSG tissue, the following signatures were significantly enriched in biopsy-positive pSS (group-IV) compared with either biopsy-negative pSS (group-III) or biopsy-negative non-SS sicca (group-I) patients: IFN-a signaling, IL-12/IL-18 signaling, CD3/CD28 T-cell activation, CD40 signaling in B-cells, double negative type-2 (DN2) B-cells, and FcRL4 + B-cells ( Figure 3A). Within group-IV, signature scores varied considerably amongst individual patients. We also assessed other signatures, such as CD40 signaling in DCs or monocytes. Only for CD40 signaling in DCs, significantly higher scores were observed in labial glands (trend for parotid glands) of biopsy-positive pSS patients (group-IV vs. group-I; data not shown). We also observed enrichment of the TLR7 signature in group-IV, with a pattern similar to the IFN-a signaling signature (data not shown), probably attributed to the large overlap in genes that make up these two signatures. None of the signatures were significantly enriched in biopsy-negative pSS compared with non-SS sicca patients ( Figure 3A).

SG Gene Signatures Reflected in Blood
We also tested gene signature enrichment in matched PBMCs, to investigate whether transcriptional changes were shared between SGs and blood. Three PBMC samples were excluded after sequencing because of a high number of reads aligning to intronic regions, indicating genomic DNA contamination. When comparing biopsy-positive pSS with non-SS sicca patients, we found that the IFN-a signaling and DN2 B-cell signatures were also enriched in pSS patients' PBMCs (Supplementary Figure 3). Upregulation of other gene signatures depends apparently on the tissue microenvironment. Importantly, IFN-a signature scores correlated significantly between blood and tissue (R 2 = 0.64 and 0.70 for PSG and LSG, respectively). Not only the PBMC IFN-a signature score, but also blood levels of MxA protein correlated significantly with IFN-a signature scores in SG tissues ( Figure 3B), indicating that blood MxA level is a biomarker for SG type-I IFN activity. SG IFN-a scores also significantly correlated with serum levels of two IFN-induced chemokines, i.e., CXCL10 and CXCL11 (Supplementary Figure 4), although not as strong as with blood MxA levels. The finding that blood MxA levels accurately reflect SG type-I IFN activity was also supported by the strong correlation between blood MxA and salivary gland MX1 transcript expression levels (R 2 = 0.66, p<0.0001 for parotid and R 2 = 0.71, p<0.0001 for labial glands). In addition to the IFN-a signature, we investigated other gene signatures elevated in SG tissue for their potential association with related serum cytokine levels. We found positive correlations between the SG DN2 B-cell score and serum CXCL13 levels, and for the CD3/CD28 T-cell activation score with serum CCL21 and sPD-1 levels (Supplementary Figure 4).

SG Transcriptome and Autoantibody Positivity
Unsupervised hierarchical clustering of pSS patients using the top-50 upregulated genes (group-IV vs. group-I) in either PSG or LSG tissue showed two main clusters: Cluster 1 with weak-tomoderate upregulation of top-50 DEGs and Cluster 2 with strong upregulation of these DEGs (Figure 4). Nearly all PSG-or LSGdefined Cluster 2 patients were anti-SSA and rheumatoid factor positive ( Figures 4A, B). Patients with FS≥3.0 in either gland were all contained within this cluster. One biopsy-positive pSS  Figure 4B). In this patient, 91% of the LSG biopsy consisted of CD45 + cells. This patient had type-III cryoglobulinemia and MALT-lymphoma in the parotid gland. Schirmer's test, UWSF and SWSF rates were all zero, reflecting the severe damage to the glands of this patient.
Lastly, we explored whether transcriptomic clusters correlated with clinical disease parameters. No significant differences in ESSDAI, ESSPRI, ESSPRI dryness, UWSF or SWSF was observed between the two clusters nor between autoantibody-positive and -negative pSS patients.

DISCUSSION
Transcriptomic analysis of paired major and minor SG biopsies of pSS patients revealed a high degree of overlap in immune pathway activity. Only tissue samples from biopsy-positive pSS were distinct from non-SS sicca patients based on the SG transcriptome. This distinction was mainly driven by infiltrating immune cells. Between biopsy-negative pSS and non-SS sicca patients, we expected to find evidence for distinct mechanisms of (intrinsic) epithelial cell dysregulation, yet the SG transcriptome of both glands was strikingly similar between these groups. For biopsy-positive pSS patients we observed that immune pathways significantly upregulated in both glands were mostly related to CD4 + T-cells and B-cells. The IFN-a gene signature was also significantly enriched in these patients and a strong correlation between the IFN-a score in PBMCs and SGs was seen, indicating that SG type-I IFN activity is reflected systemically. Our results further show that blood MxA level is a surrogate biomarker of type-I IFN activity in SGs.
We observed a strong correlation between average gene expression in PSGs and LSGs. The number of DEGs and transcript expression fold changes between biopsy-positive pSS and non-SS sicca patients were, however, higher for PSGs. This may be explained by the higher amount of infiltration in LSG vs. PSG tissue of non-SS sicca patients. Notably, in pSS patients the LSG biopsy was more often classified as positive than PSG. LSGs may be more likely to harbor non-specific infiltration, as lymphocytic foci are present in approximately 15% of healthy individuals (30,31), while for PSGs this is estimated at 5% (20). Despite these differences, Metacore pathway analysis showed that all top-20 enriched (auto)immune pathways were shared between the two different SG tissues. Importantly, immunopathways and DEGs identified in our study correspond largely to those presented in a recent transcriptomic study of LSG tissue in pSS (17).
Although many DEGs were identified in biopsy-positive pSS compared to biopsy-negative non-SS sicca patients, we found that biopsy-negative pSS patients did not have an altered gene expression profile in their PSGs and LSGs, compared to non-SS sicca patients. While not directly comparing PSGs with LSGs, we noted that genes with poorly correlated expression levels between glands were all related to SG biology, not immunopathology. Although there are indications that intrinsic epithelial cell abnormalities contribute to SG disease in pSS (32), our findings do not support this hypothesis, at least not at the transcriptomic level. It might be that our approach does not capture subtle differences between epithelial cells of biopsynegative pSS and non-SS sicca patients, or that non-SS sicca patients also exhibit epithelial abnormalities. Lastly, we cannot exclude that some non-SS sicca patients are in an early, preclinical stage of pSS, although the progression rate from non-SS sicca to pSS is estimated at only 10% over 2-3 years (33).
For biopsy-positive pSS, our study confirms the relevance of multiple SG immunopathways that have been previously published using other experimental approaches (34)(35)(36)(37)(38), although previous approaches were often focused on a single pathway and unable to measure activity of multiple pathways at the individual patient level. In the current study, upregulated gene signatures in both PSGs and LSGs of biopsy-positive pSS patients concerned IFN-a (and TLR7) signaling, IL-12/IL-18 signaling, CD3/CD28 T-cell activation, CD40 signaling in B-cells, DN2 B-cells, and FcRL4 + B-cells, with heterogeneity in signature scores amongst these patients. The signatures reveal that particularly B-and T-cells are activated, and these cells may exert pathogenic functions in the SGs. A type-I IFN signature is frequently observed in pSS patients, both in blood and SGs (26,34,39). This signature is a consequence of type-I IFN production in SGs, most likely by epithelial cells and plasmacytoid dendritic cells (37,40). IFN-a production is induced in response to (viral) RNA sensing mainly through TLR7/8 and cytosolic sensors (13). We did not find evidence for the presence of an IFN-a signature in SGs of biopsy-negative pSS patients, suggesting that this signature is tightly correlated with the amount of infiltrate, in line with literature (34). Another study showed that SS patients with FS<1.0 were less likely to have elevated expression of IFN-regulated genes in PBMCs than anti-SSA-positive SS patients with FS≥1.0 (41). We further show that there is a strong correlation between IFN-a signature scores in paired major and minor SGs and PBMCs samples. Ongoing production of type-I and other types of IFNs could maintain a pro-inflammatory state by activating immune cells and stimulating cytokine production (e.g., CXCL10, BAFF) by epithelial cells, thereby negatively affecting SG homeostasis (37,40,42,43). IFNs may also play a role in glandular T-cell activation by stimulating IL-7 production in epithelial cells (44). In our study, Tcell activation in SG lesions was reflected by an increase in the CD3/ CD28 T-cell activation score and enrichment of T-cell related pathways, with Tfh-cell dysfunction at the top. Enrichment of Tfh-cells in pSS-derived SG tissue was recently demonstrated and associated with the presence of ectopic lymphoid structures (35,45). Another gene signature enriched in SG tissue was IL-12/IL-18 signaling. IL-12, but not IL-18, was upregulated in both PSG and LSG tissue and could drive Th1-cell differentiation and type-II IFN production via STAT4 phosphorylation (46). Polymorphisms in both IL12A and STAT4 genes have been associated with pSS and patients carrying the IL12A risk allele had increased IL-12p70 serum levels (29,47).
While multiple immune cells are involved in the inflammatory response, the top-20 upregulated genes in PSG and LSG tissue of biopsy-positive pSS were mostly related to B-cells, underpinning the central role of B-cell hyperactivity in pSS pathogenesis. Combined SG microarray datasets also have shown that key driver genes in SS patients with high-grade SG inflammation were mostly related to BCR signaling and B-cell activation (19). By using SG gene signature analysis, we showed that gene signatures of activated B-cell subsets were enriched, including DN2 B-cells and FcRL4 + B-cells. Upon TLR7 stimulation, DN2 B-cells participate in extrafollicular responses in autoimmune disease patients (23). Presence of ductal epithelium-associated FcRL4 + B-cells in SGs of pSS patients was previously demonstrated (48). These cells are proliferating and phenotypically related to DN2 Bcells (24,48). We also observed increased scores for the CD40 signaling in B-cells signature in biopsy-positive pSS, in line with a recent RNA sequencing study of LSG B-cells (37). Ligation of CD40 by CD40L-expressing activated T-cells is essential for T-cell dependent B-cell activation. In pSS, chronic B-cell activation in SGs may exacerbate disease and eventually result in development of B-cell lymphoma (49).
Our results suggest that immunotherapy directed at T-cell dependent B-cell hyperactivity is particularly beneficial for biopsy-positive pSS patients. In previous studies, a positive SG biopsy has been associated with autoantibody positivity, hypergammaglobulinemia, corneal and conjunctival damage, and reduced UWSF rates (7). Furthermore, severe SG inflammation (FS≥3) has been associated with higher ESSDAI scores and higher risk of developing B-cell lymphoma (50). A positive biopsy correlates, however, poorly with symptoms of dry mouth or eyes (7), corroborating with our results that ESSPRI dryness scores were similar between transcriptomic clusters. We neither observed differences in other clinical parameters (e.g., ESSDAI, UWSF) between these clusters, emphasizing the complex relationship between local inflammation and clinical parameters assessed cross-sectionally. An apparent discrepancy between biological activity and patient-reported dryness symptoms in pSS has also been noted in several clinical trials (51)(52)(53)(54), and may be partly attributed to symptom adaptation. Trial design is further complicated by clinical and biological patient heterogeneity (55). Our results underscore our limited understanding of SG hypofunction in absence of infiltrates and question whether biopsy-negative and biopsy-positive pSS patients should receive similar treatment. For biopsy-positive pSS, our study shows enrichment of distinct gene signatures that can be therapeutically targeted, although heterogeneous expression patterns existed amongst individuals. Therefore, tissue transcriptomics may form a valuable tool for establishing precision drug therapy in pSS.

DATA AVAILABILITY STATEMENT
The salivary gland RNA sequencing datasets presented in this study can be found in online repositories. The names of the repository/ repositories and accession number(s) can be found below: https:// www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE173808.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Medical Research Ethics Committee of the University Medical Center Groningen (METc 2013.066). The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
GV, LG, SP, AV, HB, and FK contributed to conception and design of the study. AV and HB recruited the patients. GV, EH, BV, SC, VP, and FS contributed to patient material collection and/or data acquisition. LG and SH performed data quality assessment, data processing and the statistical analysis. GV, LG, SP, SM, JC, LM, and FK were involved in data analysis and interpretation of results. GV wrote the first draft of the manuscript.
LG wrote sections of the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This study was funded by Bristol-Myers Squibb (BMS). The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article or the decision to submit it for publication. The work of GV, SP, HB, and FK is supported by a Dutch Arthritis Society (ReumaNL) Long Term Project Grant (LLP-29).