Whole Exome Sequencing of Patients from Multicase Families with Systemic Lupus Erythematosus Identifies Multiple Rare Variants

In an effort to identify rare alleles associated with SLE, we have performed whole exome sequencing of the most distantly related affected individuals from two large Icelandic multicase SLE families followed by Ta targeted genotyping of additional relatives. We identified multiple rare likely pathogenic variants in nineteen genes co-segregating with the disease through multiple generations. Gene co-expression and protein-protein interaction analysis identified a network of highly connected genes comprising several loci previously implicated in autoimmune diseases. These genes were significantly enriched for immune system development, lymphocyte activation, DNA repair, and V(D)J gene recombination GO-categories. Furthermore, we found evidence of aggregate association and enrichment of rare variants at the FAM71E1/EMC10 locus in an independent set of 4,254 European SLE-cases and 4,349 controls. Our study presents evidence supporting that multiple rare likely pathogenic variants, in newly identified genes involved in known disease pathogenic pathways, segregate with SLE at the familial and population level.


Results
Exome Sequencing and Variant Filtering. We successfully sequenced the exome of the most distantly related patients with SLE (n = 5) from two well-studied multi-case SLE families from Iceland ( Fig. 1 and Table S1) 31 to a mean read depth of 33.3X across targeted coding regions (Table S2). We chose a filtering strategy favouring the best-quality single nucleotide variants (SNVs) shared by all the affected sequenced individuals from each family, with a very low allele frequency (MAF < 0.01) in an internal control population (n = 642), and with easily recognized functional consequences for protein coding genes (gain or loss of a stop codon, nonsense, missense, and essential splice sites). A total of nineteen variants passed all quality-, annotation-and frequency-based filtering criteria, and were validated by SEQUENOM and/or Sanger sequencing (Table 1). All the 19 identified variants were non-synonymous SNVs (nsSNVs) present in heterozygous state in the exome-sequenced patients. No variant with a gene truncating effect (stop gain, stop loss, nonsense, or essential splice variant) was identified. Segregation Analysis. We then performed a segregation analysis by looking at the variants co-segregating with the disease status in other affected members within each family for whom DNA was available. This analysis included five cases with SLE, one case with rheumatoid arthritis (RA) and one with multiple sclerosis (MS) in "family 6" (n = 8), as well as seven cases with SLE in "family 8" (Fig. S1).
All of the 19 identified variants were present in at least three affected members of each family including the exome-sequenced patients. Variants in KRTAP4-9 (p.D18V), KIR2DS4 (MIM 604955) (p.I255L), and SLC5A9 (p.G103R) were the most frequent among patients of both families (12, 7 and 6 patients, respectively). However, the former two variants were very frequent (MAF > 5%) among a set of 83 controls from Iceland and Sweden (referred to as Nordic controls) (Tables 1 and S3). Due to the high frequency in the Nordic controls and discrepancies across public databases, these variants were excluded from further detailed analysis.
We observed the best co-segregation with disease status in family 6 for a variant in XRCC6BP1 (also known as ATP23) (p.A229V), which was carried by six out of eight affected members including four SLE patients, one RA and one MS patient. In other words, the XRCC6BP1 variant segregated with all autoimmune diseases present in this family. The variants in SLC5A9 (p.G103R) and CHD3 (MIM 602120) (p.A1523T) were also carried by four SLE patients; the latter also by one RA patient. In family 8, five out of seven SLE relatives carried the variants in DCLRE1C (MIM 605988) (p.H283N), NOTCH1 (MIM 190198) (p.D932N), FBXL14 (MIM 609081) (p.N102H) and CLC (MIM 153310) (p.N65K) ( Table 1). Hence, there was no single variant carried by all SLE patients within each family, but rather a combination of a few rare and low frequency non-synonymous variants segregating with the disease status in most of the affected members. A large fraction of these loci was already annotated as immune-relevant, taking into account known pathways and Gene Ontology: NOTCH1, KIR2DS4, NUP214, CLC, DCLRE1C and FBXL14. Moreover, 13 of 19 variants-corresponding genes are expressed in whole blood, spleen or relevant cell lines (GTEX RNAseq data, Table 1 and Fig. 2).
Aggregate Association and Enrichment Analysis. We then tested whether any of the candidate genes identified by exome sequencing had statistical evidence of association with SLE in the general European population due to the combined effect of all rare variation within each gene (MAF < 1%) ( Table 2). For this, we used a large and independent imputed genome-wide association scan of 4,254 SLE patients and 4,349 controls with European ancestry 25 (Table 2). Each gene was analyzed using two procedures: the sequence kernel association test (SKAT) 32 and an aggregated case-control enrichment test. It is important to note that we performed targeted gene-based tests, that is, we did not test for association of rare variants neither at a genome-wide level nor tested individual variants. To be as stringent as possible, the 10 first principal components (PC's), accounting for all the significant variability due to population stratification, and genomic control (GC) were used to correct for stratification in both procedures (Figs S2 and S3). To further eliminate the potential effect of linkage disequilibrium (LD) on the computation of empirical corrected P values, tests were run using only unlinked markers by applying a very stringent LD threshold of r² < 0.1. It could be objected that such stringent filters could mask true association signals, but our rationale was that if the signals were maintained after strict correction they would strongly support a 'true positive' effect.
The exome-identified gene FAM71E1 and the adjacent EMC10 gene (MIM 614545), which is in high LD, showed significant association with both procedures, case-control enrichment and SKAT, after application of r² < 0.1 threshold and having applied the appropriate corrections for stratification and multiple testing ( Table 2). The DCLRE1C gene showed suggestive evidence of enrichment but did not remain significant after correction for multiple testing. Step-by-step scheme shows initial variant detection, validation and filtering followed by in silico functional analysis of candidate loci and analysis of candidate genes in independent population-based data.  3 and Table S3). In addition, four genes (DCLRE1C, NOTCH1, NUP214, and FAT4) were connected to specific phenotypes with immunological features through OMIM (Online Mendelian Inheritance in Man) and the Human Genome Mutation Database (HGMD) 35 . None of the nineteen variants had been reported to ClinVar 36 (https://www.ncbi.nlm.nih.gov/clinvar/) or HGMD (http://www.hgmd.org) databases (as to Jun 29th, 2017). Variants in DCLRE1C (p.H283N) and NUP214 (p.I765V) were classified as likely benign by standard algorithms, whereas the SNPDryad method predicted a likely deleterious effect. In addition, these two variants are very rare (<0.5%) in European populations (ESP6500, ExAC and 1000 genomes) and in our internal sequencing controls, which we considered as supporting evidence of a likely gene-disrupting effect (Tables 1 and S3). In contrast, the SLC5A9 nsSNV (p.G103R) (rs61746559) segregated with SLE in both families, was not common in our dataset and Nordic controls and was predicted to be damaging by several of the algorithms (Table S3), but it shows very high allele frequency in non-European populations (20% in Asians, 26% in Southern Han Chinese −1000 Genomes), which is against a gene-disrupting effect. Likewise, the variants in KRTAP4-9 (p.D18V) and KIR2DS4 (p.I255L) were rare (MAF ≤ 1%) in most European populations, but very common in the set of 85  Table 1. Segregation analysis and population frequency for exome-sequence variants. Table 1 shows the results from the segregation analysis for the variants identified by WES as the number of family members who were variant carriers affected with SLE (family 8) or any autoimmune disease (family 6), over the total number of cases in each family. Z-score is indicated for the variants located in regions linked to SLE in Iceland according to ref. 27 . MAF (minor allelic frequency) of the variants in 642 internal whole-genome sequenced European and 83 Nordic controls. The maximum MAF corresponds to the highest frequency of allele in ExAc, 1000 Genomes and GO-ESP data for global population. Variant numbers correspond to those seen in Fig. 3. Functional annotation is shown according to GeneCards SuperPathways and Gene Ontology top biological process categories. The presence of the particular gene transcripts in GTEX blood/immune samples is marked according to Fig. 3. Additional annotation of the variants shown in Table S3.
ScIEnTIfIc REPORtS | (2018) 8:8775 | DOI:10.1038/s41598-018-26274-y genotyped Nordic controls (≥5%). These common nsSNVs could be polymorphisms specific of the Asian and Nordic population, respectively. We cannot exclude that any of these are functional and/or disease variants, or that the high MAF is due to a founder effect in specific populations, but due to the high frequency and described discrepancies they were regarded as polymorphisms, unlikely to be gene-disrupting variants (Table 1, Fig. 3, and Table S3).
Gene Expression Analysis. We interrogated the potential biological relevance of the genes carrying the identified exome variants by analysing their pattern of gene expression, regulation and gene networks. Since the association of FAM71E1 extends to the nearby gene EMC10 (Table 2), and they are overlapping and transcribed in opposite directions, we included both genes in the analysis. A heat-map representing the tissue-specific RNAseq data (GTEx Project) for all genes is shown on Fig. 2 and summarized in Table S4. Most of the genes (13 out 19, or 14 out 19, taking into account expression of FAM71E1-neighbouring EMC10) were expressed in immune-relevant samples. Additionally, several genes are differentially expressed in infections and after specific pathogen exposure (Table S5). We observed no evidence of expression of KRTAP4-9, coding for hair keratin-associated protein 4-9, in any of the tissue samples analysed. The expression of PDHA2, PABPC3, and FAM71E1 was detected only in testis. EMC10 was expressed in spleen and lymphoblastoid cell lines, albeit weakly. FAT4 cadherin as well as the intestinal-specific SLC5A9 gene, although expressed in several tissues, were almost not detectable in any of the immune-relevant tissues (Fig. 2). In summary, among the exome identified genes the fraction of loci showing immune-relevant expression was high and included the EMC10 gene from the bidirectional FAM71E1/EMC10 locus. The tissue-specific expression patterns were independently corroborated using TSS activity data from the FANTOM5 project (Fig. S4).
Gene Co-expression and PPI networks. We next constructed gene networks based on the detection of common partners between the exome identified genes including EMC10. Both gene co-expression and pairwise protein-protein interactions (PPI) were taken into consideration in the networks. We first performed a family-specific analysis. The genes segregating in family 8 were part of a highly connected gene network (Fig. S5). The family-6 gene set was small; nevertheless, the genes CHD3, NUP214 and FAM8A1 were connected via at least one partner each (Fig. 5). Further on, the gene sets of both families were analysed together (Fig. 4). The resulting network included twelve genes, of which the SLE associated genes EMC10 and DCLRE1C, together with NUP214, CHD3, NOTCH1, FAM8A1, MPHOSPH8, TPRA1 and CLC showed high connectivity to each other and to other genes in the network via several partners. One interaction partner connected genes ANKRD50 and FBXL14 each, separately, to the network. Moreover, we observed strong inter-family gene connections: the family 6-specific genes, CHD3 and NUP214, were highly connected to others, including family 8-specific NOTCH1, CLC, EMC10, TPRA1, and MPHOSPH8. We independently validated the gene-gene relations by using the tissue-specific GIANT resource based on the unsupervised analysis of public RNAseq data 37 . Of note, XRCC6BP1, which was absent in the co-expression and PPI network, was connected to NOTCH1, whereas the CLC connections were not stable in the GIANT networks (Figs 5 and S6). The difference between results could probably be due to the different calculation algorithms and a  Table S3); and (right) a maximum minor allelic frequency of 1% in European populations (internal sequencing and genotyping controls, 1000 Genomes EUR, ESP6500 EurAm and ExAc Eur non-Finish populations) (See Tables S3 and 1). Genes with a significant aggregate association and enrichment of rare variants (top) have further genetic evidence of being implicated in SLE. Genes associated after removing linked variants (r 2 < 0.1), adjustment by 10 principal components and genomic control, and multiple test correction are highlighted with a star (See Table 2).  Table 2. Rare-variant Association Analysis. Gene-based case-control association analysis of unlinked rare variants (MAF < 1% and r² < 0.1) of the genes identified by exome sequencing in an independent imputed genome-wide association scan from a set of 4,254 SLE patients and 4,349 controls with European ancestry 25 . Two procedures were used, an enrichment case-control association test and the sequence kernel association test (SKAT). Correction for multiple testing was run through a bootstrapping permutation process (Pmult). All tests were corrected for stratification by adjusting for the first 10 principal components (PC) and Genomic Inflation Control (λ GC ).

Enrichment test SKAT
ScIEnTIfIc REPORtS | (2018) 8:8775 | DOI:10.1038/s41598-018-26274-y contribution of microarray-RNAseq differences in co-expression analysis. However, most of the highly connected genes coincide in the two independent approaches. Of interest, the networks contain several interaction partners previously genetically implicated in SLE, autoimmune disorders, and other immune-related disorders and traits (Table 3 and Table S6). Importantly, EMC10 has been consistently associated to the autoimmune disease primary biliary cirrhosis [38][39][40] .
Functional Enrichment. To further understand the functional relevance of the identified genes, we performed an enrichment analysis for categories using two independent tools: GeneTrail2 41 and TOPPGENE 42 . We interrogated the list of genes with nsSNVs identified in the multicase families first, and then we included the stable interaction partners identified in gene networks. We found that functional categories related to DNA metabolism and repair were significantly enriched in the exome gene list (XRCC6BP1, DCLRE1C, p (FDR-corrected) = 1.08 × 10 −4 , according to TOPPGENE), with the most significant results for V(D)J recombination (DCLRE1C, p (FDR-corrected) = 5.8 × 10 −6 ). GIANT network for tissue-specific gene connectivity. The connections are shown for "all tissues", lymph node, blood and for separate blood cell populations and correspond to one or more common partner genes between each of the studied genes. The connecting lines coloured according to the specificity. . Global pairwise gene interaction network of the exome identified genes. Gene co-expression and direct protein-protein interactions are shown as a combined network. Family 8-specific genes are shown in red ovals, family 6 -in blue ovals. The names of the genes with significant enrichment or aggregated association of rare variants in SLE cases are shown in red italic. Genes associated after removing linked variants (r 2 < 0.1), adjustment by 10 principal components and genomic control, and multiple test correction are highlighted with a star. Note the high connectivity within the network of the SLE-associated genes. Blue-filled ovals correspond to the genes with known immunity-related Mendelian disorders (OMIM data), green-filled -to the genes with published genome-wide significant associations with autoimmune and autoinflammatory disorders, yellowfilled -with other genome-wide significant immunity-related traits (Table S4).
ScIEnTIfIc REPORtS | (2018) 8:8775 | DOI:10.1038/s41598-018-26274-y Several immunity-relevant biological processes were also among the top enriched functional categories, including immune system development (NOTCH1, DCLRE1C, CLC, p (FDR-corrected) = 1.7 × 10 −6 ), and lymphocyte activation (DCLRE1C, CLC, p (FDR-corrected) = 6.5 × 10 −6 ) ( Fig. 6 and Table S7). Of note, when including mouse phenotype data in the analysis, we observed highly significant enrichment of a phenotype of absent immature B cells (DCLRE1C, p (FDR-corrected) = 5.3 × 10 −6 ). GO cellular component ontology analysis also revealed overrepresented functional groups: genes associated with DNA repair complex and nuclear chromosome. Figure 6 shows a semantic grouping of the most significant GO categories overrepresented in the list of our SLE candidate genes and their interaction partners and summarizes graphically the top enriched categories (shown also in Table S7).
Variants Putatively Linked to Gene Regulation. We observed that some of the exome variants were located in gene regulatory regions and thereby could play a role in the regulation of their respective genes. The EMC10 gene is located head-to head and shares a short bidirectional promoter with FAM71E1. According to RefSeq and GENCODE v.19 annotation, the FAM71E1 5′-UTR overlaps with the EMC10 transcription start site and the first exon. These genes are in a genomic region with very strong LD. As mentioned before, unlike FAM71E1, whose expression is restricted to testis, EMC10 has a broad expression pattern. Using Haploreg 43 for 1000 Genomes European data, we found that the FAM71E1 variant (p.L7F) (rs185418641) and its proxies are located in an open DNAse-hypersensitive chromatin site in various blood cell types. FAM71E1 is not expressed in those cell types, in contrast to EMC10.  The SLC5A9 variant (p.G103R) (rs61746559) corresponds to a immune-relevant enhancer detected in two independent functional genomics projects (FANTOM5, based of CAGE sequencing and Roadmap Epigenomics, based on histone marks). In contrast to the SLC5A9 gene, this enhancer was active in the blood, mostly in neutrophils, but not in gastrointestinal samples, similar to immune-relevant expression of candidate SLE genes.

Discussion
The pursuit of the genetic variants that explain why some individuals and their relatives have an increased risk to develop SLE is a challenging task. In our study, we focused on the identification of rare coding variants with a likely gene-disrupting effect segregating with SLE in two large pedigrees from Iceland for which we have clinical and linkage analysis data. The structure of these families might suggest a Mendelian (monogenic) pattern of inheritance, however, we found no single gene-disrupting variant segregating with the disease status in all family Figure 6. GO categories enrichment for genes with SLE-associated nsSNVs in multicase families and their interaction partners according to TOPPGENE. The overrepresented GO terms were summarized by REVIGO to representative subsets of terms using a simple clustering algorithm that relies on semantic similarity measures. The scatterplots are shown for Biological Processes and Cellular Component, the circles correspond to the cluster representatives (i.e. terms remaining after the redundancy reduction) in a two dimensional space derived by applying multidimensional scaling to a matrix of the GO terms' semantic similarities with an allowed similarity threshold = 0.9. The colour code and the size of circles reflect p-value of the overrepresented term (Table S7).
We thus looked to interrogate whether there was supporting genetic evidence for any of the family genes identified at the population level. We did so by using highly stringent rare-variant association analysis with data from a recent GWAS, a strategy previously used by us that provided support for rare variants in the RNASEH2 genes 45 . More recently, genetic association of SLE with a common variant in RNASEH2C was identified and replicated using this data in a European-Chinese meta-analysis, a variant that represented an eQTL 46 . Initial analysis based on both aggregated association and enrichment analysis of rare and low frequency variants in the identified loci supported the implication of the FAM71E1/EMC10 locus and possibly DCLRE1C. Further in silico analysis suggests a role of the variant identified by exome sequencing on the regulation of EMC10 expression. EMC10 (ER membrane complex subunit 10) codes for a protein involved in endoplasmic reticulum (ER)-associated degradation and lipid transport, but any role in SLE is completely unknown. The association analysis performed here might have been too stringent, but it allowed us to conclude with high confidence that the enrichment of rare variants at the FAM71E1/EMC10 locus is robust and not due to the effect of population stratification or LD. Regarding the other genes not associated in the GWAS data, they still might play a role in these families but not at the European population level, taking also into account possible differences between the populations studied. However, we must keep in mind that to date there are no standard methods to detect association of rare variants in complex diseases and this is a field under development with still with many challenges to address 47,48 .
Bioinformatics analysis indicated that a part of the identified variants may be predicted as having a deleterious effect, and for some variants methods did not agree. For example, the variant in DCLRE1C (DNA cross-link repair 1C) coding for the nuclear protein Artemis, was predicted as probably benign by standard algorithms but highly deleterious using a one-to-one orthologue specific approach 33 . Further analysis showed though that the detected His-> Asn change is located in a well conserved loop nearby the catalytic center and predicted DNA binding sites (Fig. S7). In addition, analysis of two variants disregarded later due to their high allelic frequency (MAF > 1%) nevertheless indicated them as worth further study: one in the neutrophil-specific enhancer in the intronic area of SLC5A9 gene and the other in the natural killer-specific KIR2DS4 gene, one of the KIR family of genes involved in the inhibition and activation of NK cell function that interact with class I MHC molecules 49 . This highlight some of the limitations of this study. First, in the absence of functional evidence, results from variant prediction algorithms can only be taken as suggestive evidence. Second, here we did not evaluate the role of common variants nor have we analysed other types of genetic variation such as structural or copy number variants.
We carried out the analysis of gene networks based on gene co-expression and protein-protein interactions and identified a functional overlap between genes, even for those segregating independently in each family. The constructed gene networks revealed among the exome genes, partners previously implicated in autoimmune diseases, including the SLE-associated genes IKZF1, IKZF3, and ITGAX 25,50,51 (Tables 3 and S6), further supporting a role for the genes identified in the disease pathogenic pathways. The importance of multiple factors and gene interactions in the genetics of complex traits is well known, and even minor fluctuations of genes expressed in the disease-relevant cells could probably contribute to the disease susceptibility, as discussed in the recent "omnigenic" hypothesis 52 . We demonstrate that even for SLE in multicase families with apparent Mendelian inheritance the underlying mechanism is complex and involves several functionally interacting genes. In contrast to family studies of Mendelian diseases, we cannot easily address issues such as incomplete penetrance nor the possibility of phenocopies, as to do so we would need to take into account all the possible genetic and non-genetic factors contributing to complex landscape of this disease.
Functional enrichment analysis of the identified genes and their interaction partners demonstrated significant overrepresentation of immunity related terms within GO categories such as DNA repair and DNA metabolism including V(D)J recombination, a critical process in the rearrangement of the T cell and B cell receptors, double strand break (DSB) repair, cellular response to DNA damage stimuli, and chromosome organization. Also, we detected a significant enrichment of genes involved in the NOTCH1 pathway ( Fig. 6 and Table S7). The DSB repair and V(D)J recombination categories were represented by DCLRE1C and XRCC6BP1 (XRCC6 binding protein 1) genes and their interacting partners. Of note, DCLRE1C recessive mutations cause Omenn syndrome (MIM 603554) (OS), a severe combined immunodeficiency (SCID) associated with increased cellular radio sensitivity due to a defect in V(D)J recombination that leads to early arrest of both B-and T-cell maturation 53 . OS displays autoimmune-like manifestations of the skin and gastrointestinal tract. SNPs in DCLRE1C have shown suggestive evidence of association with RA and T1D 54,55 . A recent functional study demonstrated that in Artemis-deficient cells type I and type III IFN signatures are elevated due to the chronic accumulation of DNA 56 .
In addition, a widely expressed heterochromatin gene MPHOSPH8 (M-phase protein 8), whose protein binds H3K9me and promotes DNA methylation 57 , was connected to DCLRE1C in all our networks, both directly and via partners. Another gene participating in chromatin regulation and highly connected to other genes in our networks was CHD3 (Chromodomain Helicase DNA Binding Protein 3). Autoantibodies against this protein are found in a subset of patients with dermatomyositis, also an autoimmune disease 58,59 . Similarly, and possibly related through XRCC6BP1, the XRCC6 gene codes for the Ku70 helicase and V(D)J recombination repair protein, a well established lupus autoantigen 60 .
Our results contribute to the growing evidence linking SLE to DNA damage and repair mechanisms (reviewed in 61 ). Increased DNA damage and radiosensitivity have been consistently reported in cells from SLE patients 61,62 . Abnormalities in V(D)J recombination in individuals with combined immunodeficiency carrying hypomorphic RAG1 pathogenic variants show manifestations of autoimmunity 63 . Importantly, abnormalities in enzymes involved in DNA metabolism have been implicated in the type I IFN response and the development of autoimmunity 64 . For example, several other genes involved in DNA repair pathways have been previously found ScIEnTIfIc REPORtS | (2018) 8:8775 | DOI:10.1038/s41598-018-26274-y associated with SLE: TREX1 (3′ repair exonuclease 1) 65,66 and X-Ray Repair Complementing Defective Repair genes including XRCC1, XRCC3, and XRCC4 67,68 .
Finally, we detected a significant enrichment of genes involved in the NOTCH1 pathway. We did not detect association of rare variants in the NOTCH1 gene itself but previous GWAS have implicated polymorphisms in this gene with several autoinflammatory diseases (Tables 3 and S5) [69][70][71] . It is worth to note also that somatic recurrent mutations in NOTCH1 and NUP214 are found in patients with hematologic malignancies for which SLE patients have an increased risk compared to the general population [72][73][74] .
In summary, we identified novel SLE susceptibility genes using exome sequencing of distantly related patients from extended pedigrees from Iceland. Taking into consideration the genetic co-segregation of variants, similar gene expression patterns, results of nsSNV protein effect prediction and gene networks modeling, we propose joint multigene mechanisms of SLE predisposition in these families. These genes highlight a role for DNA metabolism and repair in SLE pathogenesis.

Subjects and Methods
Patients and Families. We sequenced the exome of the most distantly related patients from two previously described SLE multi-case families from Iceland for which genetic linkage data was available 31 (Fig. S1). Both families have been extensively studied and have multiple cases of SLE as well as a high frequency of other autoimmune diseases 75 . All patients fulfilled the 1997 ACR classification for SLE 76 . The National Bioethics Committee (NBC) of Iceland approved the study (Approval: 02022-V4-31) and all participants gave informed consent. All experiments were performed in accordance with relevant guidelines and regulations. The details of each exome-sequenced individual and clinical data are provided in Table S1. Two SLE patients were selected from "Family 6" and three from "Family 8" (n = 5). The linkage evidence involved regions with LOD scores of 1.5-4.5. Family 8 contributed the most to the genetic linkage signals observed in the Icelandic linkage study 31 including the HLA. Family 6 did not show linkage to the HLA region 31 , thus suggesting that the genetic contribution in both families was different.
Exome Sequencing. Three micrograms of genomic DNA purified from blood of the five selected patients were enriched for coding regions using the Agilent SureSelect ® Human All Exome Target Enrichment System (38MB and 51MB kit, protocol v1.7). Deep sequencing was performed at Uppsala University on an ABI SOLiD ™ 5500xl system (Life Technologies). Colour space read correction and alignment to the Human reference sequence library (hg18) were performed at the Centro Pfizer-Universidad de Granada-Junta de Andalucía de Genómica e Investigación Oncológica (GENYO) with SOLiD ™ Bioscope Software (v.2.1, Life Technologies), obtaining a mean read depth of 33.3X across targeted coding regions (see Supplemental data and Table S2).

Variant Calling, Annotation and Filtering.
We selected only reliably mapped reads with a mapping quality (MAPQ) value over 20. PCR duplicates were removed with Picard (v1.35). Single nucleotide variants (SNVs) and indels were called by using SAMTOOLS (v0.1.10) and then exported to pileup files. Indels were called but not included in the present study because the difficulties this posed with the very short reads obtained. Information about the population frequency of Indels in public databases was also very limited at the time we selected variants for validation. SNVs in pileup file format were annotated and filtered at the Center for Human Genome Variation, Duke University, by using the Sequence Variant Analyzer (SVA) software developed by Dongliang Ge 77 . Whole genome sequence (WGS) data from 642 Caucasian unrelated individuals without any immune-related phenotype served as internal control genomes. Annotated variants in cases were annotated and filtered using the same quality control criteria and methods as the controls. Only SNVs supported by a minimum of 6 reads were included. Variants situated within UCSC Genome browser repeat masker regions were excluded. The alternative allele was compared to the chimpanzee reference allele to ensure none of the alternative alleles represent the expected ancestral allele. Tables S8 and S9 contain the summary statistics for all the shared variants filtered by SVA. The average read depth of the filtered variants in the patients was 51.44X (range 12X-180.79X).
We then selected variants by function and MAF as follows: Only protein-coding variants either introducing or removing a stop codon (stop gained and stop lost, respectively), altering a splice acceptor or splice_donor_site (essential splice site), or introducing an amino acid change (non-synonymous coding) were included. SVA uses standard Sequence Ontology (SO) definitions (http://www.ensembl.org/info/genome/variation/predicted_data.html). We considered SNVs shared between the patients of each family, with a MAF ≤ 0.5% in the control genomes. We also considered SNVs that were absent in the controls (MAF = 0), that is, those carried only by the patients (case-only variants), as well as SNVs carried in homozygous state exclusively by the patients and present in the control population only in heterozygous state with a MAF ≤ 5% (case-only-homozygous variants). We annotated variants present within the Icelandic SLE linkage regions 31 (+/−20MB) by using BEDTools. Finally, after variant validation, we also compared the allele frequencies obtained with the allele frequencies observed by the NHLBI GO Exome Sequencing Project (ESP6500) (European American population), The Exome Aggregation Consortium (European-non Finnish population) and the 1000 genomes project (Phase 3 European population) by using ANNOVAR v2012 Oct23 78 . This strict filtering strategy favours a reduction in false positives and has been widely and successfully used for the identification of disease variants in Mendelian disorders 79,80 . For further genetic analysis, the genes where the selected variants are located were referred to as candidate genes. Segregation Analysis. The identified SNVs were genotyped on a MassARRAY System (SEQUENOM) in: affected as well as healthy members from the multi-case families for whom DNA was available (Fig. 1) in 83 matched Nordic controls (n = 36 Icelandic and n = 46 Swedish controls) to determine if any of the variants were polymorphisms (MAF > 5%) specific of the North European population, and in the exome-sequenced patients as genotype controls (n = 5). We selected only individuals with a call rate per sample ≥80%. In total, six patients with SLE, one with rheumatoid arthritis (RA), one with multiple sclerosis (MS), and 1 healthy relative were included from family 6 (n = 9). Five patients with SLE, 2 individuals fulfilling 3 of 4 SLE criteria, and 4 healthy relatives were included from family 8 (n = 11). We selected only variants with a genotyping call rate ≥90% in all genotyped individuals and for which the alternative allele identified by WES was validated in the exome-sequenced patients. As we needed complete genotypes for segregation analysis, missing genotypes for the filtered variants were completed with Sanger sequencing.

Genome-Wide Association Analysis. SNP Data.
We used previously genotyped GWAS data from 5,478 individuals of European ancestry including 4,254 SLE patients and 1,224 controls genotyped as described in 25 using the Illumina© HumanOmni1_Quad_v1-0_B chip and 3,125 out-of-study controls of European origin obtained from three studies available through dbGaP with informed consent, namely the DCEG Dataset (phs000396.v1.p1; 1175 individuals), the GENIE UK-ROI Diabetic Nephropathy GWAS (phs000389.v1.p1; 903 individuals) and the High Density SNP Association Analysis of Melanoma (phs000187.v1.p1; 1047 individuals). The final data set used for aggregate/association analysis consisted of 4,212 cases and 4,065 controls (see Supplemental data).
Imputation. For each disease candidate gene, a region of interest was extended in 500,000 additional base pairs upstream and downstream, respectively, as it is known that large buffers may improve accuracy for low-frequency variants during imputation 81 . Markers within each extended region were extracted from the GWAS data for imputation with IMPUTE2 82 using the 1000 Genomes Project as reference panel. Specifically, we used 1000 Genomes Phase 3 (b37) as these haplotypes have lower genotype discordance and improved imputation performance into downstream GWAS samples, especially at low frequency variants 83 . Prior to imputation, each GWAS gene extended region was phased with SHAPEIT 81 using the EUR subpopulations as reference. A restrictive QC-filter was applied on the imputed genotypes (SNP genotyping rate ≥99%, sample genotyping rate ≥95%) without restriction of allele frequencies, in order to include both rare and low frequency variants. To ensure a highly reliable imputation, a conservative IMPUTE info_value threshold of ≥0.7 and a concordance value threshold of ≥95% for each marker were applied. We have further addressed potential bias introduced by imputation by using different association methods and keeping a very stringent significance threshold in our analyses.
Gene case-control association analysis of rare variation. Since a minor allele frequency (MAF) of 1% or more is the conventional definition of polymorphism 84 , thus we considered a MAF < 1% as 'rare variation' . We tested whether any of the genes identified by exome sequencing had statistical evidence of association with SLE in the general European population due to the combined effect of all rare variation within each gene (MAF < 1%). For this, and because there was no availability of DNA from large enough sets of patients to be sequenced, we took an alternative approach. We used a large and independent imputed genome-wide association scan from a set of 4,254 SLE patients and 4,349 controls with European ancestry 25 . Each gene was analyzed using two procedures: the sequence kernel association (SKAT) test 32 and an aggregated case-control enrichment test where adjusting a logistic regression model with a 'transformed' genetic variable equals to the sum of minor frequency alleles, below the MAF threshold, for the j markers in the gene, in each individual.
Correcting for stratification in rare variant association analysis. To be as stringent as possible, the 10 first principal components (PC) accounting for all the significant variance due to population stratification (Figs S2 and S3) and genomic control (GC) were used to correct for stratification in both procedures. The genomic inflation factor (λ GC ) was equal to 1.11 and 1.24 for 'enrichment case-control' and SKAT 10 PC-corrected tests respectively (Fig. S8). These λ GC values were used as an additional correction of the resulting inflation (GC-correction = S tatistic10PC' c_corrected/λ GC ). Without 10 PC-correction the λ GC was equal to 1.44 and 2.97 for 'enrichment case-control' and for SKAT respectively (Fig. S8). Thus, the 10 PC-correction reduced the inflation by 33% in the enrichment test, and by 174% in the SKAT test.
Correcting for multiple testing in gene case-control association analysis of rare variation. Regarding the correction for multiple testing in association of rare variants, a genomic association threshold of 10 −6 is commonly accepted (equivalent to Bonferroni correction for 19,000 to 20,000 protein encoding genes in the genome). It is also accepted that Bonferroni, although mathematically right, would be very penalizing for biological data, therefore, we opted for techniques based on permutation procedures. Our multi-test correction procedure brings together all the markers of the analysed genes into a single table; for each gene a number of markers equal to that of the analysed gene were randomly extracted from the table and its association test calculated; by repeating the procedure for N times, an empirical corrected P value (Pmult) was calculated for the analysed gene. However, in association tests that simultaneously include several markers, co-linearity due to LD between markers could potentially inflate the significance of the P value. Also, elimination of the LD by the random extraction of markers in the permutation procedure could affect the computation of empirical corrected P values. A simple correction to this is to run the association and the multiple test correction tests only with unlinked markers by applying a very stringent LD threshold of r² < 0.1. These corrected P values would depend on the observed P values used as thresholds. We have verified by linear correlation (P corrected multi-testing ~ P observed) that R² was equal to 0.99 for the enrichment tests and 0.97 for the SKAT tests. Thus the LD would not affect the correction, as expected, given the applied r² threshold.
It could be objected that applying such strict filers could mask true association signals, but given the issues related to the association of rare variants 85 , we reasoned that under all these stringent criteria any significant association signal would strongly support a real or 'true positive' association effect.  34 was used for the annotation of the potential deleterious effects of the exome variants based on the following algorithms: SIFT, PolyPhen2, FATHMM, LRT, MetaLR, MutationAccessor, MutationTester, and Provean. None of the applied methods was adapted for distinguishing effects in the paralogous proteins. Further, functional prediction of SNVs on protein-coding ENSEMBL transcripts was performed using the SNPDryad method 33 (http://snps.ccbr.utoronto.ca:8080/SNPdryad/), for which only one-to-one orthologous proteins were used for scoring. An SNPDryad score below 0.5 was assigned as possibly neutral, from 0.5 to 0.7 as possibly deleterious, and more than 0.7 as deleterious. To deal with multiple annotations, gene transcripts were first scanned for isoform-specific expression level in the GTEx Portal 86 (http://gtexportal.org). Only the main protein-coding isoforms were selected for the annotation of putative gene-damaging effects. The analysis was based on the Human Genome annotation GRC38_p3.
Tissue-specific Expression and Gene Regulation Analyses. The differential expression of genes including tissue-specific and from eQTL profiling was analysed using aggregated public microarray and RNAseq data available at GTEx 86 and MuTHER 91 . The TSS activity and gene enhancers were studied based on the FANTOM5 project CAGE data 92,93 , which contains about 900 tissue and cell specimens (http://fantom.gsc.riken.jp/5/). The regulatory enhancer elements and their tissue-specific activity were recovered from Roadmap Epigenomics data 94 (http://www.roadmapepigenomics.org/). The differential expression data were extracted from the EMBL-EBI Expression Atlas 95 , (https://www.ebi.ac.uk/gxa/about.html) Only relevant differential expression results, involving immunity, inflammation, immune-or inflammatory-stimuli and cell activation, infections, cell exposure to bacteria or bacterial components. The adjusted p-value threshold was selected as <0.05 and the absolute value of log2 fold-change as > = 2.5.
Gene and Protein Network Analyses. Gene-gene pairwise networks were constructed using two main data sources: gene co-expression data and protein-protein interactions. Gene pairs detected in two or more of the co-expression data sets were selected and included in the network analysis. Protein pairs detected in two or more sources were also included in the network analysis (See Supplemental data). Combined co-expression and protein-protein pairwise interactions were searched for direct pairs between query genes or indirect, through gene/protein partner (only one node between query genes allowed). Family-specific sub-networks were also constructed. Tissue-specific gene relations were extracted from the GIANT database 37 (http://giant.princeton.edu/; See Supplemental data).

Data availability.
In adherence with the confidentiality requirements by The National Bioethics Committee (NBC) of Iceland, individual sequences or genotype data cannot be publicly shared. All summary data generated or analyzed during this study are included in the article and its Supplementary Information files.