key: cord-0024198-kvvh7amh
authors: Teruel, María; Barturen, Guillermo; Martínez-Bueno, Manuel; Castellini-Pérez, Olivia; Barroso-Gil, Miguel; Povedano, Elena; Kerick, Martin; Català-Moll, Francesc; Makowska, Zuzanna; Buttgereit, Anne; Pers, Jacques-Olivier; Marañón, Concepción; Ballestar, Esteban; Martin, Javier; Carnero-Montoro, Elena; Alarcón-Riquelme, Marta E.
title: Integrative epigenomics in Sjögren´s syndrome reveals novel pathways and a strong interaction between the HLA, autoantibodies and the interferon signature
date: 2021-12-02
journal: Sci Rep
DOI: 10.1038/s41598-021-01324-0
sha: e685c859e8b988f9601558ef1e60d9d91be9eb3a
doc_id: 24198
cord_uid: kvvh7amh

Primary Sjögren’s syndrome (SS) is a systemic autoimmune disease characterized by lymphocytic infiltration and damage of exocrine salivary and lacrimal glands. The etiology of SS is complex with environmental triggers and genetic factors involved. By conducting an integrated multi-omics study, we confirmed a vast coordinated hypomethylation and overexpression effects in IFN-related genes, what is known as the IFN signature. Stratified and conditional analyses suggest a strong interaction between SS-associated HLA genetic variation and the presence of Anti-Ro/SSA autoantibodies in driving the IFN epigenetic signature and determining SS. We report a novel epigenetic signature characterized by increased DNA methylation levels in a large number of genes enriched in pathways such as collagen metabolism and extracellular matrix organization. We identified potential new genetic variants associated with SS that might mediate their risk by altering DNA methylation or gene expression patterns, as well as disease-interacting genetic variants that exhibit regulatory function only in the SS population. Our study sheds new light on the interaction between genetics, autoantibody profiles, DNA methylation and gene expression in SS, and contributes to elucidate the genetic architecture of gene regulation in an autoimmune population.

Differentially methylated positions associated with SS. We explored the DNA methylation patterns associated with SS comparing the whole blood DNA methylation level between 189 SS patients and 220 heathy subjects (Supplementary Table 1 ) using the Infinium MethylationEPIC BeadChip with which we could interrogate 776,284 autosomic CpG sites. In total, we observed 118 differential methylated positions (DMPs) associated with SS at a Bonferroni-corrected threshold of P < 6.4 × 10 -08 (Supplementary Table 2 ). The majority of SS-associated DMPs exhibited decreased methylation levels in SS patients compared with controls (91.5%), supporting the overall hypomethylation previously described in SS patients (Fig. 1a) 19, 20 . The 118 SS-associated DMPs were annotated to 52 unique genes and 7 intergenic regions. The majority of them fell within promoters (49.2%) and gene bodies (40.7%), and only 4.2% were located in the 3'UTR. SS-associated DMPs were mainly distributed in open sea and shore CpGs and only 5 DMPs fell in CpG islands ( Supplementary Fig. 1 ). We included an independent cohort formed by 60 SS patients and 89 healthy individuals (Supplementary Table 1 ) from whom we had DNA methylation data available from the 450 K array which contained half of the SS-associated DMPs detected in our discovery cohort. At a P < 0.05 we could successfully replicate 84.7% (Supplementary Table 2 ). The top 10 SS-associated DMPs with an average methylation difference |Δβ|> 0.1 were located within IFI44L, MX1, PARP9-DTX3L, NLRC5, IFIT1, IFIT3, IFITM1, PLSCR1, PDE7A and DDX60 genes, all of them known to be regulated by type I IFN. The most significant DMPs was the cg05696877 probe in the IFI44L gene for which an average methylation difference |Δβ|= 0.36 was observed (P = 3.5 × 10 -27 ; FDR = 2.7 × 10 -21 ) between SS and controls. Functional analyses for genes underlying SS-DMPs confirmed an enrichment of GO terms related with IFN signaling (GO:0060337; P = 3.5 × 10 -31 ; GO:0034340; P = 31.37 × 10 -30 ) as well as with defense response Table 4 ). We searched for CpGs that exhibit differences in DNA methylation variability between SS patients and controls and named these variable methylated positions (VMPs) (Supplementary Table 2 ). We observed that 80% of SS-associated DMPs also exhibited evidence of increased DNA methylation variability at a significance level of P < 0.05 in SS patients compared with healthy controls. In our replication cohort, we could replicate 83.3% of the VMPs observed. The increased DNA methylation variability associated with SS is especially pronounced for MX1 and PARP9-DXT3L genes for which we observed the largest variability differences. Increases in variabiltity might reflect epigenetic plasticity in immune-related cells and/or reflect that DNA methylation is regulated by diverse transcription factors involved in different inflammatory and immune signaling during disease progression [31] [32] [33] [34] [35] .

DNA methylation is an epigenetic mark that changes with environmental triggers as for example drugs. We investigated whether our SS-associated signatures were driven by the most common treatments in our SS patients by adjusting the linear model for treatments such as antimalarial, steroids, and immunosuppressive therapy (see details of the treatment applied to the individuals in each cohort in Supplementary Table 1 ). Our results show that most of the SS-DMPs (80.5%) remained significantly associated at our threshold P < 6.4 × 10 -08 while the reminder showed suggestive associations, of P < 1 × 10 -04 , indicating that therapy applied to the SS patients does influence the DNA methylation patterns associated with SS only in a minority of sites, with modest effects (Supplementary Table 3 ). In fact, a comparison of the effect size of the 118 significant SS-DMPs obtained in these 2 analyses, when including or not the therapy as covariate in the linear regression, showed high correlation (Pearson's correlation R = 0.99, P = 2.2 × 10 -6 ). Similar results were also observed in the replication cohort, for which we observed a high replication rate (81.0%) (Supplementary Table 3 ). Differences in DNA methylation levels and variability of CpGs located in the X-chromosome were also evaluated in females of the discovery cohort, we could not find any significant differences in females at a significance level corrected for multiple testing (data not shown). www.nature.com/scientificreports/ The IFN epigenetic signature in SS is shaped by Anti-La/SSA and associated with classII HLA genetic variation. Hypomethylation at IFN-regulated genes has been strongly associated with the presence of anti-Ro/SSA and/or anti-La/SSB autoantibodies previously 36 . Stratified analyses based on anti-La/SSB and anti-Ro/SSA autoantibodies positivity show in our data that the epigenetic signature is observed only when positive anti-Ro/SSA patients were compared with healthy subjects (Supplementary Table 3 , Fig. 2a-b) . Indeed, all associations became far from genome-wide significant (P > 0.006) when the model that includes positive and negative SS patients was adjusted by anti-Ro/SSA presence or when only patients negative for SSA were contrasted with controls (P > 0.05). Hierarchal clustering shows that DNA methylation profiles in anti-Ro/SSA negative patients resembles that observed in healthy controls rather than to anti-Ro/SSA positive patients (Fig. 2b) . While stratifying analyses based anti-La/SSB patients also yielded stronger associations in the positive group, all associations found could be explained by the presence of anti-Ro/SSA autoantibodies and nothing remained significant after adjustment for anti-Ro/SSA (Supplementary Table 3 , Fig. 2b ). Patients positive for both anti-Ro/ SSA and anti-La/SSA exhibited the strongest DNA methylation differences. See for example, in Fig. 2c , the DNA methylation difference at IFI44L gene (cg13452062) in all SS patients (|Δβ|= 0.40) goes up to |Δβ|= 0.53 when only SSA + patients are analysed, and is increased up to |Δβ|= 0.60 in the SSA + /SSB + group, and goes down to |Δβ|= 0.03 in the negative-group. This trend is conserved across all SS-associated CpGs (Fig. 2b, d, e) . Previous studies have recently shown that genetic association between HLA and SS is dependent on autoantibody profiles 36 . In order to further explore the possible link between HLA and the epigenetic IFN signature, we imputed classical HLA alleles and performed a series of logistic and linear regression models based on stratification and conditional analyses. At a significance level corrected for multiple testing we identified 2 class II HLA alleles (HLA-DQB1-0201 and HLA-DRB1-0301) associated with increased risk for SS (P < 3 × 10 -07 ) which also are associated with the Anti-La/SSA presence (P < 2 × 10 -04 ) and with the IFN epigenetic signature (P < 1 × 10 -04 ) measured as the decreased in DNAm at IFN-related genes (Table 1 ). Conditional analyses revealed that the relationship between HLA-DQB1-0201 and HLA-DRB1-0301 and the epigenetic IFN signature is lost when the presence of Anti-Ro/SSA is included in the model. We also categorized SS patients as positive or negative for the epigenetic IFN signature (see method) and observed that the association between HLA and SS is only seen in positive IFN signature SS patients (P < 8 × 10 -04 ), and not present in negative IFN signature patients (P > 0.4).

In the absence of autoantibodies, we cannot observe the HLA genetic associations with SS, and all the epigenetic signals associated with SS disappear, at least when analyzing single CpG-sites (Supplementary Table 3) . Altogether, our results suggest that there is a complex and a strong interaction between genetic variation in the HLA region, the presence of anti-Ro/SSA antibodies and the IFN-epigenetic signature observed in SS.

Differentially methylated regions associated with SS. We sought to identify differentially methylated regions (DMRs) associated with SS with the purpose of finding subtle, but consistent, methylation changes within a region that could not be detected when analyzing CpG sites by themselves. At FDR 5% threshold, we identified DMRs within 135 gene bodies, 335 promoters and 219 in CpG islands (CGI) (Supplementary Table 5 ). We found that many of the SS-associated DMPs (55.6%) were SS-associated DMRs, as for examples differentially methylated promoters for IFI44L, OAS2, RSAD2, BST2, PARP9-DTX3L, IFITM1, MX1, EPSTI1 or LY6E, and differentially methylated gene bodies of AGRN, IRF7, B2M, HERC5, ADAR and SP100, among others. Importantly, we identified 442 novel differentially methylated genes exhibiting DMRs but for which we did not find DMPs. We could not perform a strict replication of these signatures because DNA methylation was measured in a different array (450 K) in the replication sample, and this implies having lower coverage in the pre-defined regions analyzed. However, we found out that up to 45% novel DMRs, lying within 199 genes, showed robust signals in the independent sample (P < 0.05) (Supplementary Table 5 ). Among the novel genes implicated in SS we found new interferon-regulated genes, such as SAMHD1, ISG15 and XAF1; and other proteins related with the immune system such the Tumor Necrosis Factor, TNF, the TNF-receptor CD27, the chemokine receptor like protein 2, CCRL2, and the tyrosine kinase LCK. We also discovered a group of genes belonging to the HOX family, such as HOXB2, HOXD8, HOXA9, HOXA10 or HOXA4, that is implicated in transcriptional processes, as well as other transcription factors such CEBPD and GATA2. Interestingly, we discovered a large group of genes with differential methylated regions implicated in collagen metabolism such as COL11A2, COL18A1, COL27A1, COL13A1, COL23A1 and COL5A1 that were not previously identified. In addition, many DMRs, especially those that fell in CGI, were located within long non-coding RNAs (LncRNAs), as for example RP11-723C11.2, RP11-326C3.7, RP11-480D4.2 or RP11-89K21.1 (see Supplementary Table 5 ). We performed functional enrichment analyses for the set of genes showing lower DNA methylation in the SS group (64%) and separately those with increased DNA methylation (34%). On one hand, DMRs with negative signs showed an enrichment in well known pathways implicated in SS such as interferon and cytokine signaling (R-HSA-913531 and R-HSA-1280215, respectively, from Reactome database), as well as in NOD-like receptor signaling pathway (KEGG path:hsa04621), viral carcinogenesis (KEGG path:hsa05203), transcriptional misregulation in cancer (path:hsa05202), necroptosis (KEEF path:hsa04217). On the other hand, those genes showing DMRs with positive effects were found to be enriched in pathways such as Collagen biosynthesis and modifying enzymes (Reactome, R-HSA-1650814) and extracellular matrix organization (Reactome R-HSA-1474244) (Fig. 2 , Supplementary Table 6 ), which represent novel molecular pathways implicated in SS, and might indicate that the hypermethylation of these genes can be implicated in SS pathogenesis by down-regulating these molecules.

Differential expression around SS-associated epigenetic signals. In order to detect coordinated epigenetic and transcriptional changes associated with SS, we explored the possibility that the DNA methylation observed at the majority of SS-associated DMRs correlates with the expression of genes in the proximity. We identified 422 differentially expressed genes (DEGs) associated with SS when comparing RNA-seq gene Table 8 ).

Next, we investigated whether DNA methylation at DMRs and gene expression at DEGs are correlated and if this correlation is associated with SS status and could, therefore, represent coordinated effects on SS, by expression Black bar is a model that contrasted SSA-SS patients with CTRLs. Grey bar is a model that included all SS patients and was adjusted by SSA. Green bar is a model that included all SS patients and was unadjusted by SSA. Yellow line is a model that compared SSA + SS patients and CTRL. Blue line is a model that contrasted SSA + patients with CTRLs adjusted by SSB. Red bar is a model that contrasted SSA + SSB + patients and CTRLs. (c) Boxplots representing DNA methylation differences across different groups in three selected genes. R software 73 and Adobe Illustrator (https:// www. adobe. com/) was used to create figures. Table 1 . Genetic associations between HLA variation, presence of Anti-La/SSA autoantibodies and epigenetic IFN signature. β reflects the additive effect of allele dosage for different HLA alleles and P is the associated significance level. The association between HLA genetic variation, SS and SSA was determined by means of logistic regression adjusted by sex and age. The association between HLA genetic variation and epigIFN was determined by linear regression models adjusted by sex, age, cell proportions and batch effects. epigIFN refers to the epigenetic IFN signature. DNA methylation at IFI44L gene (cg13452062) was used as a proxy for epigIFN. Patients exhibiting DNAm > 0.8 were classified as negative epigIFN. Patients exhibiting DNAm < 0.8 were classified as positive epigIFN. Table 9 ). Many of the SS-eQTMs fell within the promoters of the IFN-related genes, such as IFI44L (Fig. 1c) , EPSTI1, MX1, DTX3L, PAPR9, LY6E, IFITM3, DDX60, RSAD2, PLSCR1 and ADAR. For all of them we observed that decreased DNA methylation strongly correlated with increased expression levels in SS patients, while no correlation was apparent in the healthy population. Interestingly, we observed strong coordinated effects at other non IFN-regulated genes implicated in the immune system such as CD3D, FGR, PILRA or NLRC4, and in genes with unrecognized function in SS, such as MDGA1, PM20D1 and SIRPB2 (Pearson's coefficient > 0.80). Genes showing correlated patterns of DNA methylation and gene expression were enriched mainly in immune-related gene ontologies being the most significant, immune response (GO:0006955), defense response (GO:0006952), activation of the immune response (GO:0002253) and response to stress (GO:0006950) (Supplementary Table 10 ). Our results suggest that SS-associated changes in DNA methylation can impact the transcriptional landscape in SS patients, and this could ultimately lead to alteration in the cellular function and immune response.

into the extent to which differential methylation and expression associated with SS is genetically controlled, we performed cis-meQTL and cis-eQTL analyses (see Methods). By means of linear regression models that adjust for disease conditions, we found evidence for genetic control in 52% of SS-associated DMPs and 39% of SS-associated DEGs (gene variants no farther than 1 Mb to CpG or to transcription start site-TSS). Specifically, at an FDR of 5% we found a total of 4,305 significant meQTLs that included 61 SS-DMPs and 3,508 single-nucleotide polymorphisms (SNPs). We were able to assess the replication of 1475 meQTLs that involved 31 CpGs included in the 450 K array and replicated results for 20 CpGs (64%) that were involved in 754 meQTLs (P < 0.05) (Supplementary Table 11 ). The most significant meQTLs regulate DNA methylation at the EPSTI1 gene, at an intergenic region in chromosome 3 where CCR cluster is located, and at VRK2, ADAR, IRF7 and MX1 genes. On the other hand, we identified a total of 11,399 significant eQTLs, that included 172 SS-DEGs and 10,620 SNPs, from which we could replicate 6414 eQTLs (56%) at a significance threshold of P < 0.05 and with consistent direction of effect, implicating 81 genes (Supplementary Table 12 ). Genetic variants located close to the TSS of GBP3 show the strongest association with GBP3 gene expression, followed by cis-genetic variation regulating gene expression of C3AR1 or MASTL, ETV7 and IFITM3. investigated the possibility that changes in DNA methylation and gene expression mediate genetic risk in SS, as it has been shown to occur for other autoimmune diseases 13, 33 . For that, we interrogated whether SNPs involved in meQTLs and eQTLs could be linked to the disease and showed allele frequency differences in a set of 391 SS patients and 549 healthy controls in a direction that is consistent with a mediation role of DNA methylation or gene expression (Supplementary Tables 13, 14 , Methods). At a Bonferroni-corrected significance of P < 0.0008, we detected risk variants in the HLA region of chromosome 6 that regulate DNA methylation at gene PSMB8 (Table 3 ). In this case, the minor T-allele of SNP rs7769693, at the HLA-DRB9 pseudogene, is associated with increased SS risk and with decreased methylation levels, suggesting that the variant might exert its risk by hypomethylating the PSMB8 gene. Regarding eQTL results, at a Bonferroni-corrected significant level of P < 0.0002 we detected genetic variants associated with gene expression and SS at three differentially expressed genes C2, TAP2 and PSMB9 (Table 2) . Table 2 . Genetic association of SS-meQTLs and SS-eQTLs with SS and other related SADs mediated by DNA methylation or gene expression changes. Alleles represent major allele first, and then the minor allele, which is the allele tested in each analysis. β meQTL represents the DNA methylation change with the increased in dosage of the minor allele. P meQTL corresponds to the P value from the linear regression model that regresses out the number of minor alleles for a given SNP to DNA methylation levels adjusting by age, sex, batch effects, estimated cell proportions, disease status and first genetic component. β EWAS represents the DNA methylation difference between SS and healthy controls from the epigenome-wide association study together obtained by linear regression model in which DNA methylation levels are regressed out by SS status and adjusted by age, sex, batch effects and estimated cell proportions. OR represents the Odd Ratio obtained from genetic association testing based on logistic regression modeling which the SS statuts is regressed out by number of minor alleles for a given SNP adjusted by age, sex, batch effects, estimated cell proportions and first genetic component and its corresponding. P represents the P-value obtained in the genetic associations. SAD (OR, P) represents the odd ratio and P value obtained in genetic testing for other diseases. RA = Rheumatoid Arthritis, SLE = Systemic Lupus Erythemathosus, UCTD = Undifferentiated Connective Tissue Disease, SSc = Systemic Scleroderma, PAPs = Primary anti-phospholipid syndrome. SADSnoSS = All SADs patients excluding SS. Genomic positions are based on the hg19 human reference sequence build (GRCh37). * eQTL reported in GTEx project (https:// www. gtexp ortal. org/ home/), in the case of eQTL the same SNP-gene is reported. # SNP associated with related disease phenotype in GWAS catalog (https:// www. ebi. ac. uk/ gwas/) or Open Target Genetics Portal (https:// genet ics. opent argets. org/). www.nature.com/scientificreports/ To give statistical robustness and further support to the functional and genetic link between meQTL and eQTLs with SS and with autoimmune processes, we first interrogated whether or not the novel identified SS-risk variants that reached a suggestive significance level of P < 0.05, show evidence of association with other SADs for which we had access to genotypic data: systemic lupus erythematosus (SLE), rheumatoid arthritis (RA), systemic sclerosis (SSc), mixed connective tissue disease (MCTD), undifferentiated connective tissue disease (UCTD) and primary antiphospholipid syndrome (PAPs) ( Table 2, Supplementary Table 15 ). Moreover, we searched for additional evidence of their functional role in public databases of eQTL and GWAS studies. Importantly, genetic variants regulating DNA methylation at the intergenic region of the CCR cluster in chr3 (Fig. 4a) and DNA methylation at the downstream region of IFI44 gene (Fig. 4b) show convincing evidence of their link with other SADS and have been described as eQTL for the same genes (Supplementary Table 16 ). Likewise, genetic variants regulating gene expression of the genes TRIM27 (Fig. 4c) , BTN2A2, UNC119B, GBP5 (Fig. 4d) and SLFN5 show significant associations with other SADS and a link with gene expression has been previously observed in eQTL studies for these genes (Supplementary Table 16 ). Altogether, these results support the scenario in which differences in DNA methylation and gene expression are intermediates of genetic risk for SS. This study provides a framework for the discovery of new risk loci associated with SS via alteration of regulatory landscapes.

Disease-dependent genetic effects on SS-associated differential methylation and gene expression. Finally, we hypothesized that there might exist genetic variants whose effect on molecular phenotypes depends on the specific environment of altered immune activation originated during the disease process. In order to identify such disease-dependent genetic effects or disease-interacting QTLs, we performed a gene-by-environment interaction meQTLs and eQTLs analyses in which we looked for gene variants that Table 3 . Most significant meQTLs and eQTLs exhibiting disease-dependent genetic effects. β INT represents the interaction effect between a given SNP and SS status in DNA methylation level. P INT represents the P-value obtained for the β INT in a linear regression model that adjusts for SNP, SS status, age, sex, batch effects, estimated cell proportions and the first principal genetic component. β SS.meQTL represents the DNA methylation change with the increased in dosage of the minor allele in SS population for a given SNP. P SS.meQTL corresponds to the P value from the linear regression model that regresses out the number of minor alleles for a given SNP to DNA methylation levels adjusting by age, sex, batch effects, estimated cell proportions, disease status and first genetic component in SS population. β CTRL.meQTL represents the DNA methylation change with the increased in dosage of the minor allele in the healthy control population for a given SNP. P CTRL.meQTL corresponds to the P value from the linear regression model that regresses out the number of minor alleles for a given SNP to DNA methylation levels adjusting by age, sex, batch effects, estimated cell proportions, disease status and first genetic component in the healthy control population. NA represents effects that are non significant (P > 0.05). www.nature.com/scientificreports/ interact with disease status on shaping DNA methylation and gene expression levels at SS-associated DMPs and DEGs. We considered relevant and robust disease interacting QTLs those that accomplish the following strict conditions: i) an interaction effect that passes a significance threshold of P INTER < 0.005 in the discovery cohort and a significance threshold of P INTER < 0.05 in the replication cohort, with consistent direction of the effect, ii) genetic variant associated with DNA methylation or gene expression only in the SS group (P SS.QTL < 0.05) and without evidence of genetic associations in the healthy population (P CTRL.QTL > 0.05), iii) genetic variants that show a minor allele frequency higher than 0.10 in the SS group. For DNA methylation, we observed convincing evidence for disease-interacting meQTLs at 8 SS-associated DMPs within genes TAP1, LY6E, PSMB8, STAT1, SGK269 and ATP10A (Table 3 and Supplementary Table 17 ). The most significant gene-environment interaction (ß INTER = -0.018, P INTER = 2.0 × 10 -04 ) that fell outside the HLA region involved the genetic variant rs13273708 at chromosome 8 and DNA methylation at LY6E (cg14392283) (Fig. 4a, Table 3 ). For this variant, the minor C-allele is associated with a decrease in LY6E-DNA methylation only in SS patients, but this association is not observed in the healthy population (Table 4) . A similar scenario is observed for the STAT1 gene for which we found the second most significant disease-dependent meQTLs effect (ß INTER = −0.029, P INTER = 5.47 × 10 -04 ) outside the HLA region (Fig. 5b) . Another convincing example implicates the interaction between SS status and rs1079396 in shaping DNA methylation levels at SGK269 gene (ß INTER = −0.024, P INTER = 0.003, Fig. 5c) , and the interaction with rs7169481 in shaping DNA methylation at ATP10A (cg10734665 , ß INTER = 0.022, P INTER = 0.003, Fig. 5d ). For gene expression, some IFN-inducible genes such as MX2, IFITM1, EPST1, MX1 and LY6E-DT show convincing evidence that their expression is genetically regulated in a disease-specific manner (Table 3 and Supplementary Table 18 ). The most significant diseasedependent effect was observed for MX1 gene expression (ß INTER = −0.300, P INTER = 1.8 × 10 -0 , Fig. 5e ), which is overexpressed in SS patients, and only genetically regulated within the disease population. In this case, the minor T-allele of rs9305702 is associated with decreased MX1 gene expression in SS patients, without evidence of association in the healthy population (Fig. 5e) . Other genes such as IFITM1 (Fig. 5f) , CCL2 (Fig. 4g) , NUB1 (Fig. 5h) , PLSCR1, GRIN3A and GCNT1 also show disease-specific eQTLs that regulate their SS-associated differential expression (Table 3) .

For all those disease-interacting meQTLs and eQTLs, the same trend is observed in the replication sample (Supplementary Table 17 and 18) and we found no evidence for them to be genetic regulators of DNA methylation and/or gene expression in a model that corrects for disease status and includes the whole population (Supplementary Table 17 and 18) , or in other studies that have interrogated healthy populations (Supplementary  Table 19 ) 8,10,37 . Our findings reveal a differential genetic architecture of gene regulation between SS patients and the healthy population, and suggest that the specific immunological and pathological conditions in autoimmunity 

Here, we present a comprehensive integrative large-scale analysis that serves to discovery new loci and pathways involved in SS, to recognize the importance of hypermethylation events in SS pathogenesis, and to unravel the genetic architecture of blood gene regulation in a systemic autoimmune disease. We could confirm the activation of the IFN system in SS 38 by recognizing blood hypomethylation and overexpression of a large number of genes involved in IFN signaling or specially of type I IFN-inducible genes. Type I IFNs are key immune mediators involve in viral response and in activation of immune responses 39 . Previous large-scale studies have extensively characterized the interferon signature in SS with large scale gene expression data and recently, with DNA methylation as well, in different tissues and cell types 19, 20, 22, 23, [40] [41] [42] [43] . Our stratified analyses reveal that the IFN epigenetic signature is restricted to only those SS patients that exhibit Anti-Ro/SSA autoantibodies positivity and enhanced in those that also exhibit Anti-La/SSB. It has been reported that patients positive for these autoantibodies 44 and that exhibit the IFN signature 45 are at higher risk for worse diagnosis some disease manifestations such as hypergammaglobulinemia and ongoing lymphomas. Several authors have recognized the important contribution to clinical management that stratifying patients based on gene expression based biomarkers of IFN signature 38 . Recently, a study has shown than an epigenetic-based IFN signature could be even a more suitable biomarker for patient classification 45 .

The tight relationship between autoantibodies and IFN signaling has been reported by multiple populationbased studies and experimental works in in blood, salivary glands and other target cells 22, 36, 40 . However, the mechanisms behind this close relationship are still unclear. Whether the IFN signaling have an effect on autoantibodies production or vice versa, is unresolved. On one hand, IFN signaling stimulated upon viral infection www.nature.com/scientificreports/ could generate antibodies that cross-react with autoantigens, such as Anti-Ro/SSA and Anti-La/SSB, by molecular mimicry 46 . On the other hand, it is known that the deposition of immuno-complexes formed by autoantibodies, RNA-binding proteins and material releases from apoptotic cells stimulate type I IFN production and subsequently activate IFN-regulated genes 47 . Nevertheless, autoantibodies production seems not to be enough to trigger dysregulating IFN-signaling as this is not observed in healthy subjects that exhibit autoantibodies 48 . In our study we could investigate further the genetic drivers of this strong relationship. Genetics studies have revealed that in SS there is a strong relationship between autoantibody production and class II HLA variation. Indeed, some studies have clearly show that class II HLA is only a risk factor for patients that are positive for Anti-La and/or Anti-Ro 49, 50 .In a recent SADs molecular stratification 51 , we recognized a strong relationship between class II HLA variation and a molecular cluster characterized by the IFN signature enriched in Anti-La and Anti-Ro positive patients. Beyond this previously gathered knowledge, here we have discovered that, specifically, the SS-associated class II HLA-DRB1-0301 and HLA-DQB1-0201 alleles are strongly associated with the IFN signature and with Anti-Ro/SSA and Anti-La/SSB production. Indeed, our stratified analyses revealed that class II HLA is only associated with SS and with autoantibody production in those patients exhibiting the IFN signature. Our findings point towards a complex interaction between class II HLA variation, IFN signaling and autoantibody production Other autoimmune diseases share a common IFN signature with SjS and common HLA risk variants, such as SLE and MCTD. Whether the same HLA alleles and autoantibodies are important for the IFN signature observed in other related autoimmune diseases remains to be explored. More studies that couple observational records with experimental data are needed to disentangle mechanistic routes behind these relationships and recognize potential new drug targets. Importantly, we also identified a high number of novel hypermethylation events in genes not related to the IFN pathways when applying a statistical method based on gene set enrichments 52 , a trend that was not observed when analyzing individual CpG sites. This method permits to detect associations of smaller magnitude but with consistent epigenetic patterns along pre-defined regions. Functional analyses in the group of genes exhibiting increased methylation in SS patients revealed enrichment in important pathways such as those related to the metabolism of collagen and/or implicated in extracellular matrix organization that have not been previously detected, but could explain the increased degradation of extracellular matrix structures and the significant loss of collagen observed in the lacrimal gland and other tissues in SS, and, therefore, have a key role in its pathogenesis 53 . We anticipate that future efforts in biomarker discovery will successfully recognize that these hypermethylated signals can also be of good utility for SS stratification. .

In this study, we discovered new loci associated with SS whose functional mechanisms could be the alteration of epigenetic states or gene expression profiles. The association that we have found between genetics, DNA methylation and SS in the CCR gene cluster in chromosome 3 is very interesting and represents a novel genetic risk variant. Our study shows that the minor C-allele of rs9838739, located in an active regulatory region overlapping many transcription factor binding sites and upstream the CCR gene cluster, increases risk for SS and several related systemic autoimmune diseases (SLE, MTD and UCTD) and decreases DNA methylation levels at this dense regulatory region. Previous eQTL studies from GTEx project identified that C-allele is associated with lower CCR5 expression in whole blood 37 . Interestingly, the lack of CCR5 on dendritic cells of a NOD mouse, an experimental model for SS disease and diabetes, promotes a proinflammatory environment in submandibular glands, a target and affected tissue in SS patients 54 . Furthermore, CCR5 expression is decreased on circulating monocytes from SS patients and is correlated with increased levels of inflammatory chemokines 55 . Our results, together with previous findings, suggest that genetically determined reduced methylation at CCR cluster-methylation, could contribute to SS presumably via enhanced gene expression of inflammatory chemokines. Interestingly, the same region has been strongly associated with COVID19 susceptibility and severity in several international efforts 56, 57 indicating that common molecular players may be involved in COVID19 and autoimmunity susceptibility. Other than the well described regulatory role of IRF5 risk variants 58 , we also identified that genetic variants at TRIM27, TRIM26, RLP4, NEXN, LAP3, GBP5, RBM43 and SLNF5 genes could be implicated in SS through changes in gene expression. For example, we identified a novel risk variant that we are others have demonstrated that downregulates TRIM27 37, 59 . TRIM27 is a molecule that inhibits the innate immune response and has been recently shown to negatively regulate NOD2 mediated signaling by physical interaction and degradation 60 . These data suggest that genetically determined reduced TRIM27 expression in SS patients can lead to enhanced innate responses via abnormally enhanced NOD2 activity.

Importantly, the findings of this work corroborate our hypothesis that genetic variants may only manifest genetic regulatory effects in the specific context of altered immune activity exhibited in SS patients. We identified strong meQTL effects only among SS patients and that we did not observed and have not been previously described in the general population for STAT1, LY6E and ATP10A genes 11 . STAT1 is a signal transducer and transcriptional activator that is activated in response to several cytokines as IFN-alpha, IFN-gamma or IL-6. Stimulation of monocytes and B-cells in SS patients leads to increased sensitivity of immune cells from SS patients to STAT1-activating signals that might partly explain the IFN signature observed in SS 61 In line with these results, our transcriptional data shows overexpression of STAT1 in SS patients. On the other hand, LY6E, an IFN-regulated gene that encodes for the Lymphocyte Antigen 6E, is hypomethylated and overexpressed in SS patients. Interestingly, a previous study discovered trans-eQTL effects on LY6E gene expression that is dependent on immune activation 62 , supporting that its genetic regulation is context-specific. An earlier GWAS study on cytokine responses found that genetic variants at ATP10A are associated with IFN-gamma production in response to vaccinia virus in subjects who had received the smallpox vaccine 63 . These findings, and our results, support a scenario in which inter-individual genetic variation at this gene impacts IFN-regulated gene expression only upon immune activation and via alteration of epigenetic states.

Likewise, we also identified a group of IFN-related genes which expression is genetically regulated in a disease-specific manner, such as MX1, MX2, EPST1 and IFITM1. We also discovered context-specific eQTLS Scientific Reports | (2021) 11:23292 | https://doi.org/10.1038/s41598-021-01324-0 www.nature.com/scientificreports/ for other genes with immune-related function such as CCL2 and PSLCR1. CCL2 encodes for a chemokine ligand that is secreted during inflammatory processes whose serum circulating levels are increased in SS 64 . Intriguingly, the expression of CCL2 is genetically regulated upon in vitro stimulation but it shows no genetic regulation in naïve conditions 65 . PSLCR1, the Phospholipid scramblase 1 gene, has been found to be differentially methylated and expressed in SS, SLE and MCTD 22, 36, 66 . It is a DNA-binding transcriptional activator that amplifies the IFN-mediated antiviral response and induces cytokine expression and cell proliferation. Interestingly, a previous study has shown that the genetic regulation of PSLCR1 expression is stronger upon IFN stimulation 62 . We also found disease-interacting genetic variants whose implication in systemic autoimmunity has not been well described. For example, NUB1 is an IFN-inducible gene that downregulates NEDD8, an ubiquitin-like protein.

NUB1 overexpression is known to induce inhibition of cell growth in cancer cells 67 . Its implication in autoimmune processes is uncertain but it could be related to lymphoma development. Our study has some important limitations. First of all, we focused this work on the identification of molecular changes from whole-blood derived samples, which is a mixture of different immune-cell types, many of them directly involved in the autoimmune and inflammatory processes occurring in SS. The statistical approaches that we have applied have allowed us to identify epigenetic and transcriptional signatures that are common or shared across different blood cell types. This is, on one side, an advantage to future evaluation of these signals as blood biomarkers for disease diagnosis. However, our work has not revealed molecular changes occurring in specific cell types. Previous works have recognized cell-specific epigenetic signals on data from sorted cells but on a limited sample size 43, 68 . We also could not evaluate specific epigenetic changes related to the development of some important disease phenotypes for SS, as for example lymphomas. Defective DNA methylation in SS patients with lymphomas and other cancers has been previously reported 69 and, recently a relationship between lymphoma and the DNA-methylation based interferon signature has been described 45 . However we did not have enough samples in our dataset to assess the epigenome-wide contribution in SS. Future studies that target specific disease phenotypes, that incorporate higher sample sizes, and add other cell-specific layers of molecular information will add more to the understanding of the genetic regulation and transcriptional consequences of epigenetic variation associated with SS and other systemic autoimmune diseases.

To wrap up, this study adds new pieces of evidences to the relationship between class II HLA variation, autoantibody profiles and the IFN molecular signature in SS, identifies new loci and pathways involved in SS, recognizes the importance of hypermethylation events in SS pathogenesis, and it contributes to unravel the genetic regulatory architecture of Sjögren's Syndrome, which hopefully will motivate research for the discovery of new drugs and biomarkers for SS in the near future.

Participant recruitment. Samples included in this study were recruited from the PRECISESADS study, which is a European multi center, non-randomized, and observational clinical study with recruitment performed between December 2014 and December 2018 at 19 institutions in 9 countries. PRECISESADS included patients affected by systemic autoimmune (around 400 by disease or group of diseases: rheumatoid arthritis, scleroderma or systemic sclerosis, primary Sjögren's syndrome, systemic lupus erythematosus, and primary antiphospholipid syndrome, mixed connective tissue disease and undifferentiated connective tissue disease) and 554 healthy controls 70, 71 .

A total of 558 samples were included in the core analyses of the present study (Supplementary Table 1) , that included SS patients and healthy controls for which DNA methylation and genotype data was available. All patients diagnosed with SS (N = 278) fulfilled the diagnostic criteria of the American and European community published in 2002 72 . The majority of the SS patients (94.8%) were females with a mean of age of 58.8 ± 12.4 years. Healthy individuals (N = 280), i.e., not having any history of autoimmune or infectious diseases, were included as controls, and matched to cases to the extent possible. The 74.4% of the healthy subjects were females with a mean of age of 45.1 ± 13.3 years. A large proportion of SS patients exhibited Anti-Ro/SSA positivity and/or Anti-La/SSB (70.9% and 33.7%, respectively).

Samples were divided in two sets; a discovery cohort formed by 189 SS patients and 220 healthy subjects and a replication cohort that included 60 SS patients and 89 healthy individuals. Epidemiological, clinical and serological features of the samples included in the discovery and the replication cohorts are summarized in Supplementary Table 1 . Flow cytometry data and determination of autoantibody presence was generated from serological samples as described in Barturen et al. 51 . Genome-wide DNA methylation data and differential analyses. DNA methylation data of the discovery cohort was obtained using the Infinium Methylation EPIC BeadChip (Illumina, San Diego, CA, USA) that covers more than 800,000 CpG sites. For replication cohort, we had also available genome-wide DNA meth- www.nature.com/scientificreports/ ylation data profiled obtained by Infinium Methylation 450 K BeadChip (Illumina, San Diego, CA, USA), which covers more than 400,000 CpG sites, most of them included in the EPIC array. After DNA extraction from whole peripheral blood and bisulfite conversion, the genome for each sample was amplified, fragmented and hybridized to the corresponding Illumina arrays according to the manufacturer's protocol. The quality control of samples and the normalization of the data was performed using the meffil R software 73, 74 Samples were excluded based on the detection P criteria > 99%, poor bisulfite conversion based on control dashboard check, and sex mismatches according to failed chromosome X and Y clustering. Probes were filtered out based on detection P > 0.01 in > 95% of samples. Additionally, probes located at the X and Y chromosomes were separated in different datasets to avoid gender bias. Probes with genetic variants at their CpG sites were also excluded. After applying these filtering steps we obtained 776,284 and 433,337 autosomic probes in the discovery dataset and in the replication dataset, respectively. A total of 17,530 probes located in the X chromosome passed the QC in the discovery cohort and 9282 in the replication cohort. DNA methylation was measured as a beta value ranging from 0 to 1. Zero represents an unmethylated state (0% molecules methylated at a particular sites) while 1 represents a fully methylated state (100% molecules methylated). After QC, the raw methylation beta values were background corrected and normalized using the functional normalization within the meffil R-package.

Differentially DNA methylated positions (DMPs) associated with SS were identified in our discovery cohort using a linear regression model that regresses out the SS status on DNA methylation levels at each CpG. The model was adjusted by sex, age, the first genetic principal component, the observed blood cell proportions obtained at the time of sample extraction, and batch effects as covariates. The technical variables Sentrix_ID and Sample_Plate corrected the batch effect. Data from blood cell proportions were obtained by flow cytrometry using Duraclone tubes (Beckman Coulter) by whole blood flow cytometry, as described in Jamin et al. 75 and Barturen et al. 51 (see Supplementary Note 2 for a list of investigators involved in the PRECISESADs Flow Cytometry Study Group). They were calculated as percentages of total leukocytes. Given the population stratification observed in our data, the first genetic principal component was included in the linear model as covariate. We also searched for variable methylated positions (VMP) that show DNA methylated variance differences between SS patients and healthy samples. For that, we first obtain the residuals of DNA methylation levels after correcting for all covariates in a linear regression model. Then, we searched for variance differences in DNA methylation residuals between cases and controls by applying a Levene's test that accounts for mean differences. A Bonferroni-corrected threshold of P < 6.4 × 10 -08 was established to consider DMPs and VMPs as significant. For each significant DMP the mean DNA methylation in cases and controls was calculated. The independent cohort was used to evaluate the robustness of the significant associations. The same steps were followed in the replication set. In this case, the replication significance threshold was established at P < 0.05. All statistical analyses were performed using R (v3.4.2) Software 73 . The effect of treatment and autoantibodies on the epigenetic associations found was evaluated. First, treatments and autoantibodies were, independently, included as covariates in the linear model and we compared the results of the treatment-, autoantibody-adjusted and unadjusted models. Specifically, we controlled for the use of antimalarial, steroids and immunosuppressants, since we observed that more than 10 SS patients had these drugs prescribed, and the presence of Anti-Ro/SSA and Anti-La/SSB autoantibodies. Then, stratified analyses according to the presence/absence of these treatments were performed using a linear regression model adjusted by sex, age, the first genetic principal component, the cell proportions, and batch effects. Supplementary Table 1 prevalence of treatments and autoantibodies in SS patients included in this study.

Differential DNA methylation regions (DMRs) associated with SS were estimated using mCSEA (methylated CpGs Set Enrichment Analysis) R package that implements a Gene Set Enrichment Analysis method (GSEA) to identify DMRs 52 . Briefly, the first step of mCSEA consists in ranking all the CpG probes by differential methylation and subsequently evaluates the enrichment of CpG sites belonging to the same region in the top positions of the ranked list by applying the GSEA implementation of the fgsea package. mCSEA allows to perform analysis based on promoters, gene bodies, and CGIs that are defined based on R annotation packages form Illumina for 450 K and EPIC arrays, respectively. A DMRs was considered as significant with a threshold FDR < 0.05 and with a minimum number of 5 CpGs associated to the region.

We used the mCSEA R package 52 to perform expression quantitative trait methylation (eQTM) analysis to discover significant associations between these DMRs and an expression alteration in the closer genes. In this case, the leading edge CpGs of each region is first defined and averaged for each region in each sample. Then, Pearson's correlation coefficient is calculated between each region's methylation and the proximal genes expression within 1500 base pairs upstream and downstream from the region. RNA-seq data and differential expression analysis. We analysed gene expression data available from 179 SS patients and 247 healthy subjects. RNA-seq data was obtained as 50 bp single end reads using an Illumina HiSeq 2500 sequencer and its proprietary base caller. Fastq data was aligned against human genome GRCh19 using STAR aligner 76 with GENCODE v19 [URL: https:// www. genco degen es. org/ relea ses/ 19. html] annotation as single pass alignment. Read counts were obtained using RSEM 77 . Samples were included only if they had more than 8 million reads, a RIN score > 8 and did not appear as outliers in a principal component analysis. After filtering, a total of 42,878 genes were included in the study.

Differential expression analysis was performed using a linear model that regressed out SS status on readcounts as gene expression levels adjusted by age, sex, the first genetic principal component, the cell proportions and sequencing pool (POOL) and RNA integrity number (RIN) as batch effects. A Bonferroni-adjusted P-value threshold < 1.16 × 10 -6 was established as significant to identify the differentially expressed genes (DEGs) in a discovery cohort formed by 135 SS patients and 174 healthy controls. In the replication cohort ( www.nature.com/scientificreports/ and 73 healthy controls), the same steps were followed while a replication significance threshold was established at P < 0.05. The analyses were performed using DESeq2 78 .

Genotype profiling and imputation methods. Genomic DNA from whole blood was obtained by standard methods. The samples were genotyped using HumanCore -12-v1-0-B, InfiniumCoreExome-24v1-2 and InfiniumCoreExome-24v1-3, all of them of Illumina (San Diego, CA, USA). Only genetic variants, not indels, present on all three platforms were considered for data cleaning and analysis. Quality controls (QC) were performed using PLINK v1.9 79 . Genetic markers were removed if these had a call rate < 90%, exhibited a significant differential missingness between cases and controls (P < 1 × 10 -4 ), and showed a significant deviation from Hardy-Weinberg equilibrium (P < 0.01 in controls and P < 1 × 10 -4 in cases). Variants with a minor allele frequency of less than 5% were excluded from the analysis. Samples were excluded of the study if they had a call rate < 95% and also high heterozygosity rate, i.e., they deviated 6 standard deviations from the centroid. Duplicated or related individuals were identified using identity-by-descent criteria with REAP 80 . A total of 218,947 variants passed data filtering. Samples were excluded applying a threshold of kinship coefficient < 0.25. Inference methods based on linkage disequilibrium structure was used in order to increase the number of genetic markers. Imputation was performed using the Michigan Imputation Server [URL: https:// imput ation server. sph. umich. edu/ index. html] 81 and Haplotype Reference Consortium (HRC) as reference panel [URL: http:// www. haplo type-refer ence-conso rtium. org/] 82 . We considered the imputed genotypes with a info-value (Minimac R 2 ) higher than 0.7, i.e., 70% of reliability. Imputed variants were also filtered according to the protocol described above.

In order to study the population structure and prevent population stratification the individual admixture frequency for each individual was estimated. In addition, population stratification was also analyzed by principal component analysis (PCA). Ancestry per cent for each individual was estimated using FRAPPE 83 , and a set of 2,707 independent genetic variants that maximized the differences between populations and clustering by K = 5, i. e, the 5 global populations: American, African, South Asian, East Asian and European. European published populations from 1000 Genomes phase 3 [URL: http:// www. inter natio nalge nome. org/ categ ory/ phase-3/] were included as reference panel. PCAs were also calculated using SMARTPCA from the Eigensoft software 84 and the same independent set of markers. In this analysis, only the European population from 1000 Genomes, excluding Finns, was included as reference population of our samples 85 (see supplementary Fig. 1 ). Six standard deviations from the centroid were used as a threshold to filter individuals that deviate from the main European clustering. Moreover, samples with less that 55% European ancestry were excluded from all analyses.

After genotyping QC and filtering of individuals for European ancestry classical HLA alleles were imputed for the PRECISESADS dataset (cases and controls) in the extended MHC region in chromosome 6 86 . The SNP2HLA software 87 was used for imputation using a reference panel consisting of 5,225 European individuals in the Type 1 Diabetes Genetic Consortium 88 containing data of 8,961 variants across the MHC region, and two and four digit-resolution allelic identities of the HLA class I (HLA-A, HLA-B, and HLA-C) and II genes (HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQB1, and HLA-DRB1). Genotypes were extracted as allele dosages for 298 HLA classical alleles.

Quantitative trait loci analysis and genetic associations. Genetic regulation of DNA methylation and gene expression was explored by the quantitative trait locus (QTL) approach using Matrix eQTL 89 . A linear regression model was performed to evaluate the effect of closely located SNPs (no further than 1 Mb to CpG site or TSS) with a MAF > 0.05 on DNA methylation levels of DMPs and on gene expression patterns of DEGs. In this case, the gene expression data were firstly normalized into variance stabilizing transformations (VST) using DESeq2 78 . The linear model was corrected by age, sex, disease status, the first genetic principal component, cell proportions, and the respective batch effect. A P-value threshold of false discovery rate (FDR) < 0.05 was considered as significant.

To investigate the possibility of different behavior of genetic effects between SS patients and healthy subjects in regulating DNA methylation and gene expression, meQTL and eQTL interacting analyses that included an interaction term between the considered molecular phenotype and the disease status was included in the linear regression model. Because of the lower statistical power of this analysis, we investigated only genetic variants reaching a minimum allele frequency of 10% among SS patients. We considered as significant those interaction signals that pass a suggestive P < 0.005 in the discovery cohort, that show consistent direction of effects in the replication sample and reached a P < 0.05 that could only be detected as meQTL or eQTL among SS patients (FDR < 0.05), but not in healthy subjects (P > 0.05).

Case-control genetic association analyses were conducted using PLINK v1.9 79 . Logistic regressions under the allelic additive model were performed to interrogate the association between SS diagnosis and the dosage of the minor allele in SNPs involved in meQTL and eQTLs correcting for the first genetic principal component. We used a P < 0.05 to report significant genetic associations. Moreover, we tested the genetic association of SSassociated genetic variants with other SADs comparing the allele frequencies of 261 SS, 314 SLE, 297 SSc, 332 RA, 59 PAPs, 69 MCTD and 111 UCTD patients vs. 457 healthy controls.

Logistic and linear regressions under the additive model were performed to assess the genetic association between classical imputed HLA alleles and (i) risk for SS, (ii) presence of autoantibodies, and (iii) DNA methylation levels at the IFI44L gene (cg13452062) as a proxy for the IFN epigenetic signature (epigIFN). We considered as significant those associations passing a P < 1.7 × 10 -04 that corresponds to a Bonferroni-adjusted significance considering the number of HLA alleles tested. We also categorized SS patients as those exhibiting epigenetic hypomethylated IFN signature (DNA methylation-beta value < 0.8 at cg13452062) and those being negative www.nature.com/scientificreports/ for the IFN signature (DNA methylation-beta value > = 0.8) and then performed HLA genetic associations in a stratified-fashion between the epigIFN positive and negative groups.

Bioinformatics tools. Significant DMPs were annotated to genes and gene locations according to annotation files provided by Illumina, from which we obtained a list of unique differentially methylated genes. Enrichment analyses based on gene ontology (GO) terms were performed using the web server ConsensusPathDBhuman [URL: http:// cpdb. molgen. mpg. de/ CPDB] 90 that integrates different types of functional interactions from 30 public sources in order to assemble a more complete and a less biased picture of cellular biology. Currently, ConsensusPathDB contains metabolic and signaling reactions, physical protein interactions, genetic interactions, gene regulatory interactions and drug-target interactions in human, mouse, and yeast.. We performed gene-set overexpression analyses using as a background the total list of autosomal genes covered by the Illumina EPIC array, a minimum set-size of 2 next-neighbors, and the following pathway databases: Reactome, KEGG, Wikipathways and Biocarta. For gene ontologies, we selected gene ontologies levels of 2 and 3 categories, and pathways, and both biological and molecular functions. Representations of enriched pathways were performed using cluster profiler R packages 91 .

The cohort datasets generated and analyzed during the current study are available upon request through ELIXIR platform.

Received: 13 January 2021; Accepted: 19 October 2021

Sjögren's syndrome

Primary Sjögren's syndrome

Genetics of Sjögren's syndrome in the genome-wide association era

Genome-wide association studies in Sjögren's syndrome: What do the genes tell us about disease pathogenesis?

Genetics of systemic lupus erythematosus and Sjögren's syndrome: An update

Variants at multiple loci implicated in both innate and adaptive immune responses are associated with Sjögren's syndrome

A genome-wide association study of global gene expression

Systematic identification of trans eQTLs as putative drivers of known disease associations

Finding the missing heritability of complex diseases

Disease variants alter transcription factor levels and methylation of their binding sites

Genomic and phenotypic insights from an atlas of genetic effects on DNA methylation

Integrative genomic deconvolution of rheumatoid arthritis GWAS loci into gene and cell type associations

Epigenome-wide association data implicate DNA methylation as an intermediary of genetic risk in rheumatoid arthritis

Mapping eQTLs with RNA-seq reveals novel susceptibility genes, non-coding RNAs and alternative-splicing events in systemic lupus erythematosus

DNA methylation mapping identifies gene regulatory effects in patients with systemic lupus erythematosus

Linking genetic variation with epigenetic profiles in Sjögren's syndrome

Perceptions of epigenetics

Genetic and environmental impacts on DNA methylation levels in twins

Genome-wide DNA methylation patterns in naïve CD4+T cells from patients with primary Sjögren's syndrome

Genome-wide DNA methylation analysis in multiple tissues in primary Sjögren's syndrome reveals regulatory effects at interferon-induced genes

Epigenome-wide association studies for systemic autoimmune diseases: The road behind and the road ahead

Peripheral blood gene expression profiling in Sjögren's syndrome

Gene expression profile in the salivary glands of primary Sjögren's syndrome patients before and after treatment with rituximab

Overlap between differentially methylated DNA regions in blood B lymphocytes and genetic at-risk loci in primary Sjögren' s syndrome

Expression quantitative trait locus analysis in systemic sclerosis identifies new candidate genes associated with multiple aspects of disease pathology

Cell Specific eQTL Analysis without Sorting Cells

Genetics of gene expression in primary immune cells identifies cell type-specific master regulators and roles of HLA alleles

Genomic modulators of gene expression in human neutrophils

Identification of context-dependent expression quantitative trait loci in whole blood

Genetic regulatory effects modified by immune activation contribute to autoimmune disease associations

Epigenetic and Transcriptional Variability Shape Phenotypic Plasticity

Differential variability analysis of gene expression and its application to human diseases

Epigenome-wide comparative study reveals key differences between mixed connective tissue disease and related systemic autoimmune diseases

Increased DNA methylation variability in type 1 diabetes across three immune effector cell types

Increased DNA methylation variability in rheumatoid arthritis-discordant monozygotic twins

Shared and unique patterns of DNA methylation in systemic lupus erythematosus and primary Sjögren's syndrome

Genetic effects on gene expression across human tissues

Interferons in Sjögren's syndrome: genes, mechanisms, and effects

Immunomodulatory functions of type I interferons

Activation of IFN pathways and plasmacytoid dendritic cell recruitment in target organs of primary Sjögren's syndrome

MicroRNA-146a contributes to abnormal activation of the type I interferon pathway in human lupus by targeting the key signaling proteins

Epigenetic signatures of salivary gland inflammation in Sjögren's syndrome's. Arthritis Rheumatol

Cell-specific epigenome-wide DNA methylation profile in long-term cultured minor salivary gland epithelial cells from patients with Sjögren's syndrome

How immunological profile drives clinical phenotype of primary Sjögren's syndrome at diagnosis: analysis of 10,500 patients (Sjögren Big Data Project)

DNA methylation-based interferon scores associate with sub-phenotypes in primary Sjögren's syndrome

Autoimmunity and infection in Sjögren's syndrome

Activation of the type I interferon system in primary Sjögren's syndrome: a possible etiopathogenic mechanism

Serum type I interferon activity is dependent on maternal diagnosis in anti-SSA/Ro-positive mothers of children with neonatal lupus

Genetic and clinical basis for two distinct subtypes of primary Sjögren's syndrome

In primary Sjögren's syndrome, HLA class II is associated exclusively with autoantibody production and spreading of the autoimmune response

Integrative analysis reveals a molecular stratification of systemic autoimmune diseases

MCSEA: detecting subtle differentially methylated regions

Increased degradation of extracellular matrix structures of lacrimal glands implicated in the pathogenesis of Sjögren's syndrome

Lack of CCR5 on dendritic cells promotes a proinflammatory environment in submandibular glands of the NOD mouse

Lack of CCR5 on dendritic cells promotes a proinflammatory environment in submandibular glands of the NOD mouse -Wildenberg -2008

Genetic mechanisms of critical illness in COVID-19

TITLE Mapping the human genetic architecture of COVID-19 by worldwide meta-analysis

A common haplotype of interferon regulatory factor 5 (IRF5) regulates splicing and expression and is associated with increased risk of systemic lupus erythematosus

The interferon-signature of Sjögren's syndrome: How unique biomarkers can identify underlying inflammatory and immunopathological mechanisms of specific diseases

TRIM27 Negatively Regulates NOD2 by Ubiquitination and Proteasomal Degradation

Cytokine-induced STAT1 activation is increased in patients with primary Sjögren's syndrome

Innate immune activity conditions the effect of regulatory variants upon monocyte gene expression

Genome-wide analysis of polymorphisms associated with cytokine responses in smallpox vaccine recipients

Circulating cytokines in primary Sjögren's syndrome determined by a multiplex cytokine array system

Innate immune activity conditions the effect of regulatory variants upon monocyte gene expression

An integrative cross-omics analysis of DNA methylation sites of glucose and insulin homeostasis

NUB1, an interferon-inducible protein, mediates anti-proliferative actions and apoptosis in renal cell carcinoma cells through cell-cycle regulation

Genome-wide DNA methylation patterns in naive CD4+ T cells from patients with primary Sjögren's syndrome

MTHFR gene variants and non-MALT lymphoma development in primary Sjogren's syndrome

Towards the taxonomy of human disease

Integrative analysis reveals a molecular stratification of systemic autoimmune diseases

Classification criteria for Sjögren's syndrome: a revised version of the European criteria proposed by the American-European Consensus Group

R: A language and environment for statistical computing. R Foundation for Statistical Computing

Meffil: efficient normalization and analysis of very large DNA methylation datasets

Multi-center harmonization of flow cytometers in the context of the European 'PRECISESADS' project. Autoimmun

STAR: Ultrafast universal RNA-seq aligner

RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

PLINK: a tool set for whole-genome association and population-based linkage analyses

Estimating kinship in admixed populations

Next-generation genotype imputation service and methods

A reference panel of 64,976 haplotypes for genotype imputation

Fast model-based estimation of ancestry in unrelated individuals

Principal components analysis corrects for stratification in genome-wide association studies

An integrated map of genetic variation from 1,092 human genomes

A high-resolution HLA and SNP haplotype map for disease association studies in the extended human MHC

Imputing amino acid polymorphisms in human leukocyte antigens

Overview of the MHC fine mapping data

Matrix eQTL: ultra fast eQTL analysis via large matrix operations

ConsensusPathDB: Toward a more complete picture of cell biology

clusterProfiler 4.0: A universal enrichment tool for interpreting omics data

Pôle de pathologies rhumatismales systémiques et inflammatoires

Servicio Andaluz de Salud

29 Serviço de Imunologia EX-CICAP, Centro Hospitalar e Universitário do Porto, Porto, Portugal. 30 IMIBIC, Reina Sofia Hospital

Pôle de Pathologies Rhumatismales Inflammatoires et Systémiques, Institut de Recherche Expérimentale et

39 Laboratorio di Analisi Chimico Cliniche e Microbiologia -Servizio di Citofluorimetria, Fondazione IRCCS Ca' Granda Ospedale Maggiore Policlinico di Milano

Zuzanna Makowska and Anne Buttgereit are employees of BAYER AG. All other authors have no competing interests. Reprints and permissions information is available at www.nature.com/reprints.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.