ARTICLE https://doi.org/10.1038/s41467-020-16327-0 OPEN An integrated multi-omics approach identifies the landscape of interferon-α-mediated responses of human pancreatic beta cells Maikel L. Colli 1 ✉, Mireia Ramos-Rodríguez2,3, Ernesto S. Nakayasu 4, Maria I. Alvelos1, Miguel Lopes1, Jessica L. E. Hill5, Jean-Valery Turatsinze 1, Alexandra Coomans de Brachène1, Mark A. Russell5, Helena Raurell-Vila2,3, Angela Castela1, Jonàs Juan-Mateu1, Bobbie-Jo M. Webb-Robertson4, Lars Krogvold6, 1234567890():,; Knut Dahl-Jorgensen6, Lorella Marselli 7, Piero Marchetti7, Sarah J. Richardson 5, Noel G. Morgan 5, Thomas O. Metz 4, Lorenzo Pasquali 2,3,8 & Décio L. Eizirik 1,9,10 Interferon-α (IFNα), a type I interferon, is expressed in the islets of type 1 diabetic individuals, and its expression and signaling are regulated by T1D genetic risk variants and viral infections associated with T1D. We presently characterize human beta cell responses to IFNα by combining ATAC-seq, RNA-seq and proteomics assays. The initial response to IFNα is characterized by chromatin remodeling, followed by changes in transcriptional and transla- tional regulation. IFNα induces changes in alternative splicing (AS) and first exon usage, increasing the diversity of transcripts expressed by the beta cells. This, combined with changes observed on protein modification/degradation, ER stress and MHC class I, may expand antigens presented by beta cells to the immune system. Beta cells also up-regulate the checkpoint proteins PDL1 and HLA-E that may exert a protective role against the auto- immune assault. Data mining of the present multi-omics analysis identifies two compound classes that antagonize IFNα effects on human beta cells. 1 ULB Center for Diabetes Research, Medical Faculty, Université Libre de Bruxelles, Brussels 1070, Belgium. 2 Endocrine Regulatory Genomics, Department of Experimental & Health Sciences, University Pompeu Fabra, 08003 Barcelona, Spain. 3 Endocrine Regulatory Genomics Laboratory, Germans Trias i Pujol University Hospital and Research Institute, Badalona, Spain. 4 Biological Sciences Division, Pacific Northwest National Laboratory, Richland, WA 99352, USA. 5 Institute of Biomedical & Clinical Science, University of Exeter Medical School, Exeter EX2 5DW, UK. 6 Division of Pediatric and Adolescent Medicine, Faculty of Medicine, Oslo University Hospital, Oslo, Norway. 7 Department of Clinical and Experimental Medicine, Islet Cell Laboratory, University of Pisa, 56126 Pisa, Italy. 8 Josep Carreras Leukaemia Research Institute (IJC), Badalona, Barcelona, Catalonia, Spain. 9 WELBIO, Université Libre de Bruxelles, Brussels, Belgium. 10 Indiana Biosciences Research Institute, Indianapolis, IN, USA. ✉email: mcolli@ulb.ac.be NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 T ype 1 diabetes (T1D) is a chronic autoimmune disease in human islets10,14 are significantly enriched in immune-related leading to pancreatic islet inflammation (insulitis) and pathways, including type I and II interferon regulation/signaling progressive beta cell loss1. Type I interferons (IFN-I), a (Supplementary Fig. 2b). Next, we performed a Rank–Rank class of cytokines involved in antiviral immune responses2, are Hypergeometric Overlap (RRHO) analysis (which estimates the involved in insulitis. Viral infections are a risk factor associated similarities between two ranked lists15) comparing the log2 fold- with T1D development3 and individuals at risk of T1D show a change (FC) ranked list from RNA-seqs of EndoC-βH1 cells and type I interferon signature4. The type I interferon, interferon-α human islets (IFNα-treated vs untreated) against an equally (IFNα), is expressed in islets of T1D patients5, and antibodies ranked list of genes obtained from RNA-seq of purified primary neutralizing different isoforms of IFNα prevent T1D development beta cells16 from T1D and healthy individuals (Supplementary in individuals with polyglandular autoimmune syndrome type 16. Fig. 2c and Supplementary Data 1). There was a significant Exposure of human pancreatic beta cells to IFNα recapitulates intersection between upregulated genes induced by IFNα in both, three key findings observed in human insulitis, namely HLA class EndoC-βH1 cells (362 overlapping genes) and human islets (850 I overexpression, endoplasmic reticulum (ER) stress and beta cell overlapping genes), and genes induced by the local pro- apoptosis7. inflammatory environment affecting primary beta cells from Combination of genome-wide association studies (GWAS)8 T1D individuals (Supplementary Fig. 2d, f). We also compared and studies using the ImmunoChip9 have identified around 60 these two IFNα-treated datasets against beta cells from T2D loci associated with the risk of developing T1D. Transcriptomic patients17, a condition mostly characterized by metabolic stress18. studies revealed that >70% of the T1D risk genes are expressed in By contrast with the observations made in beta cells from T1D human pancreatic beta cells10, and many of these genes regulate individuals, there was no statistically significant correlation innate immunity and type I IFN signaling11. between IFNα-regulated genes in EndoC-βH1 cells and human Type I IFN signaling is often cell-specific, an effect mediated by islets and the gene expression profile present in T2D beta cells differences in cell surface receptor expression, and activation of (Supplementary Fig. 2e, g). downstream kinases and transcription factors12. Thus, and con- sidering the potential relevance of this cytokine to the patho- genesis of T1D, it is crucial to characterize its effects on human IFNα induces early changes in chromatin accessibility. The beta cells. To define the global impact of IFNα on human beta ATAC-seq experiments demonstrated that INFα induces early cells, we presently performed an integrative multi-omics analysis changes in chromatin accessibility, with >4400 regions of gained (ATAC-seq, RNA-seq and proteomics) of IFNα-treated human open chromatin regions (OCRs) detected at 2 h, which decreased beta cells to determine the early, intermediate and late responses to 1000 regions by 24 h (Fig. 1b and Supplementary Data 2); only to the cytokine. The findings obtained indicate that IFNα pro- nine regions had loss of chromatin accessibility (Fig. 1b). Most of motes early changes in chromatin accessibility, activating distant the OCRs at 24 h were already modified at 2 h (fast response), and regulatory elements (RE) that control gene expression and pro- only 10% of OCRs were specifically gained at 24 h (late response). tein abundance. IFNα activates key transcription factors (TFs), The gained OCRs were mostly localized distally to gene tran- including IRF1, which act as a mediator of the crosstalk between scription starting sites (TSS) (Supplementary Fig. 3a) acting, beta cells and immune cells via the expression of the checkpoint therefore, as potential regulatory elements. These regions are proteins PDL1 and HLA-E. Furthermore, IFNα induces modules evolutionary conserved (Supplementary Fig. 3b), and enriched for of co-expressed mRNA and proteins that physically interact and transcription factors (TFs) binding motifs (Supplementary have relevance to T1D pathogenesis. The integration of high- Fig. 3c), including islet-specific TFs binding sequences. coverage RNA-seq and ATAC-seq indicates regulatory gene To assess whether changes in chromatin remodeling were networks and reveals that alternative splicing and different first associated with variations in gene expression, we first quantified exon usage are key mechanisms expanding the repertoire of the frequency of ATAC-seq regions gained or stable in the mRNAs and proteins expressed by stressed beta cells. Finally, proximity (40 kb window centered on the TSS) of genes with mining the modules of co-expressed genes and the IFNα beta cell differential mRNA expression (up/down/non-regulated or non- signature against the most recent catalogs of experimental and expressed) (Supplementary Data 2). There was a higher clinical drugs identifies two potentially interesting therapeutic proportion of upregulated genes associated with gained OCRs targets for future trials. in comparison to stable regions at each time point analyzed (Fig. 1c). Moreover, the number of gained OCRs was associated with changes in both the proportion (Fig. 1d) and the intensity Results (Supplementary Fig. 3e) of transcript induction (Supplementary IFNα modifies beta cell mRNA expression similarly to T1D. Fig. 3d, see Methods for more information). There was also a We performed a time course multi-omics experiment combining minor association between the number of stable regions and ATAC-seq, RNA-seq and proteomics of the human beta cell line upregulated mRNAs at 2 h (Supplementary Fig. 3e), likely due to EndoC-βH1 exposed or not to IFNα. The data were integrated to the activation of already nucleosome-depleted regions ahead of determine the dynamics of chromatin accessibility, gene/tran- cytokine exposure19. Consistently with these results, there was an script expression and protein translation, respectively (Fig. 1a). increase in the frequency of upregulated proteins coded by genes We also performed RNA-seq of 6 independent human pancreatic proximal to gained OCRs (Fig. 1e). Likewise, there was a islet preparations exposed or not to the cytokine at similar time progressive increase in IFNα-induced protein abundance depend- points (Supplementary Fig. 1a). To assess whether our in vitro ing on the number of linked gained open chromatin regions model is relevant for the in vivo islet inflammation (insulitis) in (Fig. 1f). T1D, we took two approaches: (1) Examine whether candidate There was a strong correlation between upregulated mRNAs genes for T1D expressed in human islets are involved in IFN and induced proteins (r2: 0.66 and 0.65 at, respectively, 8 and signaling (Supplementary Fig. 2a); and (2) Compare our in vitro 24 h, p < 2.2 × 10−16) (Fig. 1g, first column), but a much lower data of IFNα-treated EndoC-βH1 cells and human islets with similarity between downregulated mRNAs and proteins (Fig. 1g, available RNA-seq data of human beta cells from T1D patients. In second column). Gene ontology analysis of differentially line with previous findings suggesting a role for IFNs on the abundant proteins upon IFNα treatment identified several pathogenesis of T1D13, we found that T1D risk genes expressed biological processes involved in the pathogenesis of T1D, such 2 NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 ARTICLE a Untreated 2h 8h 24 h Experimental design 2.0 IRF1 (MA0050.1) 1.5 Bits 1.0 0.5 0.0 1 2 3 4 5 6 7 8 9 10 11 12 Data integration Cut-site probability 0.4 IFNα 0.2 0.0 –100 –50 0 0 50 100 ATACseq Dist. to motif (bp) ATACseq RNAseq + + + RNAseq EndoC-βH1 cells RNAseq Proteomics + Proteomics b c d 2h RNA type ATAC-seq –log10 p value 40 ATAC 2 h ATAC 2 h ATAC 24 h ATAC 24 h 30 ATAC 2 h ATAC 24 h stable gained stable gained 1.00 20 1.00 0.75 0.75 RNA 2 h RNA 2 h 10 0.50 0.50 0.25 0.25 0 0.00 0.00 –2 0 2 4 1.00 1.00 ATAC-seq log2FC Proportion Proportion 0.75 0.75 RNA 8 h RNA 8 h 24 h 0.50 0.50 ATAC-seq –log10 p value 40 0.25 0.25 0.00 0.00 30 1.00 1.00 RNA 24 h 0.75 0.75 RNA 24 h 20 0.50 0.50 10 0.25 0.25 0.00 0.00 0 gained stable gained stable 0 1 >1 0 1 >1 0 1 >1 0 1 >1 –2 0 2 4 ATAC-seq type # OCRs of same type ATAC-seq log2FC e f g Up-regulated (8 h) Down-regulated (8 h) 2 2 0 r = 0.66 2 r = 0.0065 –16 –16 p value < 2.2 × 10 p value < 2.2 × 10 # of gained OCRs 1 Protein type >1 4 1 IFITM2 ATAC 2 h ATAC 24 h BST2 WDR45 FBF1 HRNR SYT2 ATAC 2 h ATAC 24 h 1.00 0 ADAT3 AMIGO3 LONRF2 CERK HPDL CMTM4 2 MX1 1.00 IFIT3 IFIT2 0.75 SP110 SAMD9L UBE2L6 –1 0.75 OAS2 0.50 0 8h SP100 8h 0.50 –2 0.25 0.25 0 5 10 15 20 –3 –2 –1 0 Density Proportion 0.00 0.00 1.00 1.00 Up-regulated (24 h) Down-regulated (24 h) 2 2 0.75 0.75 r = 0.65 2 r = 0.065 –16 –16 p value < 2.2 × 10 p value < 2.2 × 10 24 h 0.50 24 h 0.50 4 1 0.25 0.25 PCDHA10 0.00 OAS2 0 NT5DC2 JMJD7 IGDCC4 0.00 2 MX1 BST2 CERK CMTM4 PADI2 PRR11 gained stable gained stable ISG15 IFIT3 SAMD9L SP100 LONRF2 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 IFIT2 SP110 UBE2L6 ATAC-seq type TRIM22 –1 Proteomics log2FC IFITM2 0 XAF1 –2 0 5 10 15 –3 –2 –1 0 Fig. 1 Exposure of EndoC-βH1 cells to interferon-α promotes changes in chromatin accessibility, which are correlated with gene transcription and translation. a EndoC-βH1 cells were exposed or not to IFNα (2000 U/ml) for the indicated time points and different high-throughput techniques were performed to study chromatin accessibility (ATAC-seq, n = 4), transcription (RNA-seq, n = 5) and translation (Proteomics, n = 4). b Volcano plot showing changes in chromatin accessibility measured by ATAC-seq. Open chromatin regions indicated as gained (red) or lost (blue) had an absolute log2 fold- change (|log2FC|) > 1, and a false discovery rate (FDR) < 0.05. The regions that did not reach such threshold were considered “stable” (gray). c, d Frequency of upregulated, downregulated or stable transcripts in the vicinity (<20 kb transcription start site (TSS) distance) of one or multiple open chromatin regions (OCRs) as classified in b. e Frequency of differentially abundant proteins in the vicinity (<20 kb TSS distance) of gained or stable open chromatin regions. f Distribution of IFNα-induced changes in protein abundance among upregulated proteins based on the number of linked gained OCRs. g Correlation between RNA-seq and proteomics of EndoC-βH1 cells exposed to INFα. The x axis represents the mRNA log2FC. The most upregulated (log2FC > 0.58, FDR < 0.05) and downregulated (log2FC < −0.58, FDR < 0.05) mRNAs are filled in red and blue, respectively. The y axis indicates the proteomics log2FC. The proteins most upregulated (log2FC > 0.58, FDR < 0.15) or downregulated (log2FC < −0.58, FDR < 0.15) are represented by red and blue borders, respectively. mRNAs and proteins not meeting these criteria were considered equal-regulated (gray fill and border, respectively). NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 a b 2h 24 h 10 IRF1 (MA0050.1) STAT1 (MA0137.1) 2.0 2.0 1.5 Bits 1.5 1.0 Bits 1.0 0.5 A 0.5 0.0 0.0 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 13 14 0.6 0.4 8 Cut-site probability 0.3 Cut-site probability 0.4 0.2 0.2 0.1 6 0.0 –100 –50 0 0 Dist. to motif (bp) 50 100 0.0 –100 –50 0 0 Dist. to motif (bp) 50 100 STAT2 (MA0517.1) 2.0 1.5 4 IRF1 (MA0050.1) Bits 1.0 0.5 0.0 2.0 B RNA-seq log2FC 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1.5 Bits 1.0 0.5 0.0 1 2 3 4 5 6 7 8 9 10 11 12 0.6 0.5 Cut-site probability 0.4 0.4 Cut-site probability 2 10 0.3 C 0.2 0.2 0.1 0.0 8 D –100 –50 0 0 Dist. to motif (bp) 50 100 0.0 –100 –50 0 0 Dist. to motif (bp) 50 100 0 2h 8h E 24 h IRF1 (MA0050.1) 6 2.0 1.5 Bits 1.0 F 0.5 0.0 1 2 3 4 5 6 7 8 9 10 11 12 4 RNA-seq log2FC –2 0.4 2 Cut-site probability 0.2 0 2h 8h 24 h –4 0.0 –100 –50 0 0 Dist. to motif (bp) 50 100 –2 –4 –6 –6 c NT 0.5 h 1 h 2h 4h 8h 24 h d siCT siIRF1 e siCT siSTAT1 siSTAT2 p < 0.0001 f siCT siIRF1 g siCT siSTAT1 siSTAT2 p < 0.0001 h siCT siSTAT1+2 1.2 p < 0.0001 1.2 p < 0.0001 1.2 1.2 1.2 p < 0.0001 p = 0.0001 p < 0.0001 p < 0.0001 p < 0.0001 p < 0.0001 p < 0.0001 PDL1/Beta-actin (FOLD) PDL1/Beta-actin (FOLD) PDL1/Beta-actin (FOLD) PDL1/Beta-actin (FOLD) PDL1/Beta-actin (FOLD) IFNα – – + – + – + – + – + – + p-STAT2 0.8 24 h 0.8 0.8 24 h 0.8 0.8 2h total STAT2 p = 0.0006 0.4 0.4 0.4 0.4 0.4 Tubulin 0.0 0.0 0.0 0.0 0.0 NT 0.5 h 1 h 2h 4h 8h 24 h NT IFNα NT IFNα NT IFNα NT IFNα NT IFNα IFNα – – + – + – + – + – + – + p-STAT1 i siCT j siCT siSTAT1 k siCT l siCT m siCT siSTAT1 total STAT1 1.2 sIRF1 p = 0.0009 siSTAT2 p = 0.037 1.2 sIRF1 p = 0.0025 siSTAT2 p = 0.039 siSTAT1+2 1.2 1.2 1.2 p < 0.0001 p < 0.0001 p = 0.0052 p = 0.007 HLA-E/Beta-actin (FOLD) HLA-E/Beta-actin (FOLD) HLA-E/Beta-actin (FOLD) HLA-E/Beta-actin (FOLD) p < 0.0001 HLA-E/Beta-actin (FOLD) p < 0.0001 Tubulin 0.8 0.8 0.8 0.8 p = 0.0012 0.8 24 h 24 h 2h NT 0.5 h 1 h 2h 4h 8h 24 h IFNα – – + – + – + – + – + – + 0.4 0.4 0.4 0.4 0.4 IRF1 Tubulin 0.0 NT 0.0 NT 0.0 NT 0.0 NT 0.0 NT IFNα IFNα IFNα IFNα IFNα Fig. 2 IRF1, STAT1 and STAT2 regulate IFNα-induced transcription and the expression of checkpoint proteins. a The regulatory paths summarize the temporal patterns of the differentially expressed genes (DEG) detected by RNA-seq (|log2FC|> 0.58 and FDR < 0.05, n = 5) (evaluated by DREM25). The x axis represents the time and the y axis the mRNA log2FC. Each path corresponds to a set of co-expressed genes. Split nodes (circles) represent a temporal event where co-expressed genes diverge in expression. In blue are the TFs upregulated at the respective time points of the RNA-seq that may regulate the pathways. b IFNα promoted TFs footprint deepening in open chromatin regions (OCR) associated to genes from the indicated DREM pathways. OCRs were associated to the nearest gene TSS with a maximum distance of 1 Mb. Previously annotated TF matrices79 were used to identify differential DNA-footprints induced by IFNα (blue lines = untreated cells, red lines = IFNα (24 h), dashed lines = reverse strand, continuous line = forward strand, Methods, n = 4). c Time course profile of STAT1, STAT2 and IRF1 protein activation in EndoC-βH1 cells exposed to IFNα (representative of four independent experiments). d–m EndoC-βH1 cells were transfected with an inactive control siRNA (siCT) or previously validated7,24 siRNAs targeting IRF1 (siIRF1), STAT1 (siSTAT1), STAT2 (siSTAT2) or STAT1 plus STAT2 (siSTAT1 + 2). After 48 h the cells were exposed to IFNα The values were normalized by the housekeeping gene β-actin (mRNA) and then by the highest value of each experiment considered as 1 (for h and m (n = 3); for e–g, i, j and l (n = 4); for d, k (n = 5)), ANOVA with Bonferroni correction for multiple comparisons (d–m). Values are mean ± SEM (d–m). Source data are provided as Source Data file. as antigen processing and presentation, responses to viruses, of differentially expressed genes (DEG) from all RNA-seq time apoptosis and NK/T-cell responses (Supplementary Fig. 4a, b); points (Supplementary Data 3) was analyzed using the dynamic groups of genes associated to protein modification and degrada- regulatory events miner (DREM) model23. This approach iden- tion were also present (Supplementary Fig. 4a, c). Furthermore, tified six patterns of co-expressed genes (Fig. 2a); 5 out of 6 genes related to endoplasmic reticulum (ER) stress, another post- pathways had an early peak of induction (2 or 8 h), which then transcriptional mechanism that downregulates translation of decreased or remained stable until 24 h (Fig. 2a). The model many mRNAs20, were also upregulated by IFNα at both the compared the frequency of TF binding sites in the gene pro- mRNA and protein levels (Supplementary Fig. 4d). These moters between divergent branches of co-expressed genes, findings are in line with our previous observations7 and were assuming that these TFs are responsible for the observed differ- confirmed here in independent samples for two key ER stress ences in gene expression profiles (Fig. 2a). This was compared markers, namely the transcription factor ATF321 and the ER with the TF occupancy determined by assaying the protection of chaperon HSPA5 (also known as BiP/GRP78)22 (Supplementary the bound sequence to ATAC-seq transposase cleavage (foot- Fig. 4e–h). ER stress often decreases translation, which may print) (Supplementary Fig. 5a and Methods). There were foot- explain the weak association observed between mRNA and prints for the transcription factors IRF1, STAT1 and STAT2, protein expression in downregulated mRNAs and proteins which were deepened upon IFNα exposure in pathway B (which (Fig. 1g). had the highest transcriptional upregulation at 2 h) and for IRF1 in two independent pathways, namely B and D at 24 h (Fig. 2b). Western blot analysis confirmed the activation of these TFs IRF1, STAT1 and STAT2 are key regulators of IFNα signaling. (Fig. 2c). STAT1 and STAT2 phosphorylation peaked between 0.5 To identify the key transcription factors involved, the expression and 1 h and then returned to near-basal levels at 24 h, while IRF1 4 NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 ARTICLE peaked later, at 4–8 h decreased by 24 h (Fig. 2c); these findings alpha cells as compared to beta cells (Fig. 3j). This may help to support the observed TF footprint profiles (Fig. 2b). There was explain why alpha cells are more resistant to the immune assault also a close correlation between DEGs induced by IFNα in RNA- in T1D. seq of EndoC-βH1 cells and in human pancreatic islets (Sup- plementary Fig. 1b; p < 2.2 ×10−22 at 2, 8 and 24 h), which mRNA and protein modules regulated by interferon-α. We resulted in a similar pattern of gene activation under the control integrated the RNA-seq and proteomics data (using all the of analogous TFs (Supplementary Fig. 1c and Supplementary samples from both 8 and 24 h) using the weighted correlation Data 4). network analysis package (WGCNA)28. The heatmaps of the Individual DREM pathways usually regulate specific biological topological overlap matrix from each dataset with module processes (GO) (Supplementary Fig. 5b, 1d). Among them, was assignment are shown in Fig. 4a. There were initially 32 eigengene the term “Regulation of immune responses” (Supplementary modules of mRNAs and 27 of proteins, which were merged Fig. 5b). This pathway comprises several genes involved in the (considering a dissimilarity threshold of 0.25) reducing the crosstalk between beta cells and the immune system, such as numbers of mRNA and protein modules to 8 and 7, respectively PDL1 (CD274), an immune checkpoint protein expressed in the (Supplementary Fig. 7a–c). The quality of these modules was islets of T1D individuals24, and a second co-inhibitory molecule, determined using a combined score of density and separability HLA-E, recently identified as potential target for cancer measures (Methods)29, which indicated that they were well- immunotherapy25 (Fig. 2d–m). defined (Zsummary > 10) (Supplementary Fig. 7d). WGCNA ana- By using a previously validated siRNA targeting IRF124, we lysis of the RNA-seq of human islets exposed to IFNα identified obtained around 60% knockdown (KD) of INFα-induced IRF1 well-defined modules of mRNAs (Supplementary Fig. 8a–d), protein and mRNA expression at 2 and 24 h (Supplementary similar to the ones identified in EndoC-βH1 cells exposed to the Fig. 5c–f). IRF1 silencing led to a significant decrease in IFNα- cytokine (Supplementary Fig. 8e). To focus on central modules induced PDL1 and HLA-E mRNA expression (Fig. 2d, f, i, k). induced by IFNα exposure, we selected only the differentially Silencing of IRF1 also decreased IFNα-induced upregulation of expressed genes (DEG) (Supplementary Data 3) and abundant the chemokines CXCL1 and CXCL10, the HLA-I component proteins (DAP) (Supplementary Data 5) in each eigengene beta-2-microglobulin (B2M) and the suppressor of cytokine module, representing 49% of the protein-coding DEGs and 89% signaling 3 (SOCS3) (Supplementary Fig. 5d, f). Small interference of the DAPs, and then examined the overlap between these RNAs targeting STAT1 (siSTAT1) or STAT2 (siSTAT2) pro- datasets. There was a significant overlap between five modules of moted >70% KD of their respective proteins and mRNAs, mRNAs and proteins (minimum of 10 elements in common, (Supplementary Fig. 5g–j). Inhibiting STAT1 or STAT2 alone FDR < 0.05) (Fig. 4b). The two main new modules, called #1 and partially blocked the induction of PDL1 and HLA-E at 2 h #2 (Fig. 4c), were composed of highly correlated mRNAs and (Fig. 2e, j), but led to a paradoxical increase in PDL1 and HLA-E proteins (Supplementary Fig. 7e, g) predominantly upregulated expression at 24 h (Fig. 2g, i), which is probably due to a by IFNα at both 8 and 24 h (Supplementary Fig. 7f, h). Module #5 compensatory increase in expression of the non-targeted STAT24. also had significantly correlated members (Supplementary In line with this, double KD of STAT1 + STAT2 led to Fig. 7i), but enriched in downregulated mRNAs/proteins at both downregulation of both PDL1 and HLA-E (Fig. 2h, m). STAT2 8 and 24 h (Supplementary Fig. 7j). Interestingly, there was sig- inhibition decreased the 2 h expression of IFNα-induced CXCL1/ nificant enrichment of ATAC-seq gained OCRs in module #2 10, SOCS1 and MX1, whereas STAT1 KD only prevented CXCL10 (Fig. 4d). They were enriched for TF binding motifs including induction (Supplementary Fig. 5h). At 24 h only 2 out of 4 genes both the pro-inflammatory motifs ISRE / IRF and the islet- remained partially inhibited by siSTAT2 (Supplementary Fig. 5j), specific transcription factor FOXA2 (Fig. 4e). whereas double KD of STAT1 + STAT2 prevented IFNα-induced To identify the gene regulatory network (GRN) of module #2, gene upregulation at 24 h in most cases (Supplementary Fig. 5k). we integrated information from two sources: (1) literature-based Exposure of FACS-purified human beta cells (Supplementary collection of TF-target interactions30, and (2) the present de novo Fig. 6a–c) to IFNα confirmed the upregulation of genes related to TF binding motifs and their predicted targets (Supplementary antigen presentation (HLA-I), antiviral responses (MX1, MDA5), Fig. 9a). This allowed us to add information from cis-regulatory ER stress (CHOP), immune cells recruitment (CXCL10) and elements (in orange) acting on the IFNα-induced GRN in human checkpoint regulators (PDL1) (Supplementary Fig. 6d). beta cells (Supplementary Fig. 9b). A similar approach was used The checkpoint protein PDL1 is overexpressed in beta cells for modules #1 and #5, but considering only data from the from people with T1D24, and we presently evaluated the literature (Supplementary Fig. 10a, c). expression of another checkpoint protein, i.e. HLA-E25. IFNα The PPI network InWeb InBio Map31 was used to assess the upregulated HLA-E mRNA expression in EndoC-βH1 cells presence of protein–protein interaction (PPI) networks in the (Fig. 3a), dispersed human islets (Fig. 3b) and FACS-purified different modules. This generated networks of interacting human beta cells (Fig. 3c) and augmented HLA-E protein proteins for modules #1, #2 and #5 (Fig. 4f and Supplementary expression in both EndoC-βH1 cells (Fig. 3d) and human islets Fig. 10b, d) and allowed the recognition of protein communities (Fig. 3e), with peak at 24 h. The inhibitory effects of HLA-E on (grouped by colors) that regulate specific and common biological immune cells require its expression on the cell surface or its functions (Fig. 4f and Supplementary Fig. 10b, d). Module #2, secretion26. Flow cytometry confirmed that IFNα increases which presents the higher number of connections, showed an surface HLA-E expression (Fig. 3f, g, Supplementary Fig. 5l), enrichment for several key biological processes activated by IFNα but there was no HLA-E release to the supernatant (Supplemen- and relevant for the pathogenesis of T1D, including cellular tary Fig. 5m). HLA-E mRNA expression was upregulated by 8- response to viruses, antigen processing and presentation via MHC fold in human islets of donors with recent-onset T1D in the class I, inflammatory and acute phase responses (Fig. 4g). DiViD study27 and HLA-E protein expression was significantly increased in insulin-containing islets, but not in insulin-deficient islets, of T1D individuals in comparison to healthy individuals Interferon-α changes the alternative splicing landscape. The (Fig. 3h, i). HLA-E expression was present in both beta and alpha present high-coverage RNA-sequencing (>200 million reads) cells (but not delta cells; Supplementary Fig. 5n) in the islets of allowed the detection of ~47,000 splicing variants, with IFNα- people with T1D, with a predominance of expression among induced 343 differentially expressed transcripts (DETs) at 2 h, NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 a 1.2 b 1.2 p < 0.0001 c 1.2 p = 0.0006 p < 0.0001 p < 0.0001 HLA-E/Beta-actin (FOLD) HLA-E/Beta-actin (FOLD) HLA-E/Beta-actin (FOLD) 0.8 0.8 0.8 p < 0.0001 p = 0.0014 0.4 0.4 0.4 0.0 0.0 0.0 2h 8h 24 h 8h 24 h 24 h d 8h 24 h e 8h 24 h f g IFNα – + – + IFNα – + – + HLA-E HLA-E Tubulin Tubulin Control IFNα Mean Fluorescence Intensity 1.2 1.2 800 300 p = 0.015 p < 0.0001 p < 0.0001 Unstained HLA-E/α-Tubulin (FOLD) HLA-E/α-Tubulin (FOLD) p < 0.0001 Control 600 IFNα 0.8 0.8 200 Count 400 0.4 0.4 100 200 0.0 0.0 0 0 8h 24 h 8h 24 h 100 101 102 103 104 105 24 h HLA-E h T1D ICI HLA-E DAPI HLA-E GLUC DAPI HLA-E GLUC INS DAPI T1D IDI HLA-E DAPI HLA-E GLUC DAPI HLA-E GLUC INS DAPI No diabetes HLA-E DAPI HLA-E GLUC DAPI HLA-E GLUC INS DAPI i 150 j ns HLA-E GLUC INS HLA-E expression MFI (AU) p < 0.0001 100 50 0 No diabetes T1D ICI T1D IDI and 1690 and 1669, respectively at 8 and 24 h, with predominance 8 h and 2344 events at 24 h (FDR < 0.05) (Fig. 5b). The most of upregulated transcripts (Fig. 5a and Supplementary Data 6 and frequent AS event modified by IFNα was cassette exons (CEx), 7). Considering all the DETs, 4%, 32% and 32% were exclusively with predominantly increased exon inclusion (represented by modified at 2, 8 and 24 h, respectively, indicating a predominance ΔPSI > 0.2, FDR < 0.05) (Fig. 5c). An example of a cassette exon of intermediary to late transcriptional changes induced by IFNα. showing increased inclusion upon IFNα treatment is the gene Next, we evaluated the frequency of each individual splicing OASL (Fig. 5d, e), an antiviral factor targeting single-stranded events (with an absolute difference in percent spliced-in (|ΔPSI|) RNA viruses such as picornaviruses32. Exposure to IFNα for 24 h > 0.2) regulated by IFNα at 8 and 24 h. There were 3140 events at increased exon 4 inclusion in both EndoC-βH1 cells and human 6 NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 ARTICLE Fig. 3 HLA-E is overexpressed in pancreatic islets of T1D individuals. EndoC-βH1 cells (a, d), human islets (b, e) or FACS-purified human beta cells (c) were exposed (gray bars) or not (black bars) to IFNα for the indicated time points and HLA-E mRNA (a–c) and protein (d, e) evaluated. The values were normalized by the housekeeping gene β-actin (mRNA) or α-tubulin (protein) and then by the highest value of each experiment considered as 1 (for a (n = 4); b (n = 3 (8 h), n = 5 (24 h)); c (n = 4); d (n = 4) and e (n = 2 (8 h), n = 4 (24 h)), ANOVA with Bonferroni correction for multiple comparisons (a–e)). f, g HLA-E cell surface expression was quantified in EndoC-βH1 cells by flow cytometry. Histograms (f) represent changes in mean fluorescence intensity (MFI). The MFI values (g) were quantified at baseline and after 24 h exposure to IFNα (n = 4, two-sided paired t-test). Values are mean ± SEM (a–g). h Immunostaining of HLA-E (green), glucagon (red) and insulin (light blue) in representative islets from individuals with or without diabetes. The top and middle panels represent an insulin-containing islet (ICI) and insulin-deficient islet (IDI) from T1D sample DiViD 3, and the lower panel represents an islet from a control donor (EADB sample 333/66). DAPI (dark blue). Scale bar 20 μm. i The MFI analysis of HLA-E expression. 30 ICIs from 6 independent individuals with T1D (5 islets per individual), 20 IDIs from 4 independent individuals with T1D (5 islets per individual), and 30 ICIs from 6 independent individuals without diabetes (5 islets per individual) were analyzed. Values are median ± interquartile range; ANOVA with Bonferroni correction for multiple comparisons, AU (arbitrary units), ns = (non-significant). j Higher magnification image demonstrating that HLA-E (green) localizes predominantly to alpha cells in a T1D donor islet (glucagon (red); insulin (light blue)) but is also expressed in beta cells, as indicated in h and j. Scale bar 30 μm. Source data are provided as Source Data file. islets (Fig. 5d). In line with this, the protein encoded by the Table 1)40 in IFNα-treated EndoC-βH1 cells as compared to the isoform OASL−001 (which retains exon 4) displayed a higher remaining proteins (Fig. 5j). IFNα-induced upregulation in comparison with the protein encoded by the isoform OASL−002, which has exon 4 exclusion (Fig. 5e). Interestingly, the isoform OASL−001 has antiviral IFNα induces increased alternative transcription start sites. activity, whereas the isoform 002 lacks the ubiquitin-like domain The usage of alternative transcription start (TSS) sites is another required for this response (Supplementary Fig. 11A)33. mechanism that generates different transcripts from the same Intron retention is an important mechanism of gene expression gene41. We used the SEASTAR pipeline42 for the computational regulation, promoting nuclear sequestration of transcripts or identification and quantitative analysis of first exon usage. This cytoplasmatic degradation via nonsense-mediated decay34. There approach recognized >250 events of alternative first exon (AFE) was a predominance for intron removal after 24 h (represented by usage occurring in 166 different genes at 8 h, and >130 events of ΔPSI < −0.2, FDR < 0.05), but not at 8 h (Fig. 5f). To understand AFE usage in 88 genes at 24 h (Fig. 6a). In agreement with this, how this impacts protein translation, we compared changes in 118 and 64 alternative promoters (±2 kb around FE TSS) detected protein abundance among three categories of ΔPSI. Genes by SEASTAR at 8 and 24 h, respectively, overlapped peaks of TSS presenting intron removal had a significant increase in protein identified by the FAMTOM5 Consortium43. Among these genes expression after IFNα exposure for 24 h in comparison to those was the 5′-nucleotidase cytosolic IIIA (NT5C3A), a negative with intron retention (ΔPSI > 0.2, FDR < 0.05) or with non- regulator of IFN-I signaling44. This gene had two AFEs identified significant intron changes (ΔPSI −0.2–0.2 or FDR > 0.05) by the SEASTAR modeling. In untreated condition (controls), (Fig. 5g). there was a higher usage of the proximal first exon (FE), present There were clear variations in the mRNA expression of several in the isoforms NT5C3A−001 and 002 in beta cells (Fig. 6b, well-known RNA-binding proteins (RBPs)35 upon IFNα exposure upper panel). After INFα exposure, however, there was increased (Fig. 5h, left panel), but the impact on the respective proteins was usage of the distal FE from the transcript NT5C3A−004 (ΔPSI: less pronounced (Fig. 5h, right panel). We focused on a group of 0.71 (8 h)/0.65 (24 h), both FDR < 0.001), which is supported by IFNα-modified RBPs at both mRNA and protein levels after 24 h, the cap analysis of gene expression (CAGE) of TSSs45 (Fig. 6b, and mapped their RNA-binding motifs among upregulated and upper panel). This was confirmed in independent samples of downregulated alternative exons. In support of a biological role for EndoC-βH1 cells and human islets using specific primers (Fig. 6b, these RBPs on alternative exon splicing, there was an enrichment lower panel). Exon Ontology analysis46 indicated that this FE of their binding motifs in regions controlling alternative cassette shift probably has functional impact, since the distal FE lacks exon inclusion/exclusion (Fig. 5i). To further study some of these both the endoplasmic reticulum (ER) retention signal and the findings, we first reproduced the IFNα-induced downregulation of transmembrane helix (Supplementary Fig. 11b), enabling its two RBPs, ELAV-like protein 1 (ELAVL1) and heterogeneous encoded protein to remain in the cytosol where NT5C3A acts44. nuclear ribonucleoprotein (HNRNPA1), by using specifics siRNAs Next, we compared the frequency of gained OCRs among (Supplementary Fig. 12a, e). Next, we evaluated whether this alternative promoters. As the SEASTAR pipeline mainly inhibition reproduced the changes induced by IFNα in the exon recognizes non-redundant FEs, we evaluated alternative promo- usage of four-and-a-half LIM domain protein 1 (FHL1) and ters identified by both the SEASTAR pipeline and the FAMTON5 Caprin Family Member 2 (CAPRIN2) (Supplementary Fig. 12b, f) database of alternative TSSs45 (Supplementary Methods). We two potential targets of, respectively, ELAVL136 and HNRNPA137. thus identified 198 and 51 gained OCRs present in alternative Silencing these RBPs promoted changes on exon usage (Supple- promoter regions at 2 and 24 h, respectively. Characterization of mentary Fig. 12c, g) that were similar to the ones observed after the IFNα-induced alternative promoters presenting a major gain IFNα treatment (Supplementary Fig. 12b, f). This is especially in chromatin accessibility pointed to the T1D risk gene RMI247. relevant in the context of the IFNα-induced exon exclusion FHL1, At gene level, there was only a ~1.4-fold upregulation of RMI2 which decreases the expression of transcripts coding for the expression, but at the transcript level there was a >60-fold protein FHL1A (Supplementary Fig. 12d), an isoform described as increase in two isoforms, RMI2−002 and −004. Visualization of a key host factor for the replication of the RNA virus the RMI2 locus combined with ATAC-seq and RNA-seq peaks Chikungunya38. indicated that the isoform RMI2−004 gained chromatin acces- RBPs can also control gene expression by blocking RNA sibility in its promoter leading then an increase in mRNA translation, as described for the Fragile X Mental Retardation 1 expression (Fig. 6c). Data from CAGE analysis45 and RNA (FMR1) gene39. Indeed, there was a significant downregulation of polymerase II ChIP-seq of another human cell type exposed to previously validated bona fide targets of FMR1 (Supplementary IFNα48 (Fig. 6c, lower part) confirms the presence of the RMI2 NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 a RNA seq Proteomics b 1 2 RNA seq clusters 5 3 –log p value 4 4 3 5 2 6 1 0 7 8 1 2 3 4 5 6 7 Proteomic clusters c d e Consensus Name HOMER matches p value Targets (%) Background (%) 200 ** Chromatin type IRF1/2/3/7/8, ISRE, 800 ISRE/IRF STAT1::STAT2, PRDM1 1 × 10–25 28.6% 0.4% both up-regulated gained Number of elements both down-regulated stable MEIS1/2/3, TGIF1/2, –14 MEIS 1 × 10 25.0% Number of regions 150 1.2% PKNOX1/2 RNA up, prot down 600 ZBTB ZBTB12 (0.66) 1 × 10–13 10.7% 0.0% 100 400 ZEB1 (0.65) 1 × 10–13 ZEB (Zf) 8.9% 0.0% FOXA FOXA2 (0.65) 1 × 10–13 26.8% 1.8% 50 200 RUNX (Runt) RUNX2, KLF1/4/9 1 × 10–12 8.9% 0.0% 0 0 FOXJ FOXJ1 (0.65) 1 × 10–12 14.3% 0.2% (R M rot4 #1 (R M rot5 #2 (R M rot1 #3 (R M ot2 4 ro #5 #1 #3 #4 #2 #5 MEF2 r # 1 × 10–12 N o ) N o ) N o ) N o ) ) MEF2D (0.69) 30.4% 3.1% , P le , P le , P le , P le , P le t1 e e e e e (MADS) ul ul A1 du A1 du A2 du A2 du A5 du ul ul ul od od od od od N Mo 1 × 10–12 M M M M M ARID ARID3A 7.1% 0.0% (R f g TAPBP PRKD2 TAP2 HLA-F LGALS8 TAP1 HLA-E LNPEP PSMB8 CTH PDIA4 SORL1 PSMB10 PSMB9 PLEKHA7 GNB4 DNAJA1 STX17 CDH1 HLA-C response to virus NAPA BAZ1A response to type I interferon CLN3 MAX EPSTI1 GSDMD MIER1 CXCL1 CNP response to interferon-alpha CASP4 RFX5 MT2A PARP9 NMI CXCL10 IRF9 regulation of I-kB kinase/NF-kB signaling DAB2IP PARP14 USP25 Size RSAD2 FBXO6 IFI35 CASP7 STAT1 STAT2 antigen processing and presentation via MHC-I RNF135 50 TNFSF10 OPTN IRF7 RNF213 DTX3L UBE2L6 TRAFD1 membrane fusion 100 FADD PML FMR1 TRIM21 IRF3 150 TRIM25 HERC5 JAK1 protein ADP-ribosylation TBK1 CASP8 RBCK1 LONRF1 FAM46A 200 TRIM5 MAVS AZI2 RIPK1 PARP10 PLSCR1 response to tumor necrosis factor TRAF3 IFIH1 DDX60 –log FDR TNFRSF11A regulation of transposition 10.0 IFIT3 DDX58 7.5 USP21 UTRN inflammatory response 5.0 EIF2AK2 2.5 DHX58 acute-phase response BST2 0.0 antigen receptor-mediated signaling MYD88 PIK3AP1 pyrimidine nucleoside catabolic process IRF1 IRF2 response to exogenous dsRNA GTPBP2 IFIT5 APOBEC3G ERAP1 MAP2 ATXN7 TDRD7 ETV6 IRF8 extrinsic apoptotic signaling pathway POLH UBA7 IFIT1 SP100 ADAR ETV7 0.1 0.2 0.3 0.4 APOBEC3F ERAP2 KATNA1 AGRN CDK17 IFIT2 MX1 gene ratio ISG15 AGO2 LAP3 HELZ2 GBP1 USP18 ZC3HAV1 OAS3 NUB1 OASL OAS2 Fig. 4 Weighted correlation network analysis (WGCNA) identifies IFNα-regulated mRNA and protein modules. a Heatmap representation of the topological overlap matrix. Rows and columns correspond to single genes/proteins, light colors represent low topological overlap, and progressively darker colors represent higher topological overlap. The corresponding gene dendrograms and initial module assignment are also displayed. b Identification of modules presenting significant overlap (FDR < 0.05 and a minimum of 10 members in common) (green border) between differentially expressed genes (DEG) and their translated differentially abundant proteins (DAP). c Composition, number of elements and type of DEG and DAP present in each of the significantly overlapping modules. d ATAC-seq-identified open chromatin regions at 2 h were linked to gene transcription start sites (TSSs) in a 40 kb window. These genes and their open chromatin regions were associated to the modules of DEG and DAP. The enrichment for gained open chromatin regions was then evaluated in each module. (** represents a p-value = 0.002343, one-sided χ2 test). e De novo HOMER motifs present in the ATAC-seq regions overlapping module #2 as described in Methods. The unadjusted p-values were obtained using the hypergeometric test from the HOMER package77. f The protein–protein interaction (PPI) network of module #2 was done using the InWeb InBio Map database31. Enriched proteins (FDR < 0.05 and minimum number of connections = 5, represented as squares) were identified and added to the network if they were not already present. Red fill identifies upregulated proteins, blue fills downregulated proteins and gray fill equal-regulated. Colored regions delimitate communities of proteins, as described in Methods. The wordcloud next to each community presents their enriched geneRIFs terms. g The biological processes (GO) overrepresented in module #2 summarize the main findings observed in IFNα-treated human beta cells. The present results were based on RNA-seq (n = 5) and proteomics (n = 4) data of EndoC-βH1 cells. 8 NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 ARTICLE a f 1500 Number of transcripts 2.5 3 2.0 density 1000 density 1.5 2 24 h 8h Down-regulated Up-regulated 1.0 1 500 0.5 0.0 0 −1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0 0 ΔPSI ΔPSI 2h 8h 24 h b g p < 0.0001 PSI neg 1477 206 573 446 438 p < 0.0001 (47%) (6.5%) (18.3%) (14.2%) (13.9%) 8h without RI Events A3SS 1118 161 409 322 334 A5SS PSI pos (47.7%) (6.9%) (17.5%) (13.7%) (14.3%) MXE RI 24 h CEx 00 05 10 15 20 25 5 .0 0. 0. 0. 0. 0. 0. –0 log2 FC c h 24 h 24 h 8h 8h ANKHD1 ANKHD1 3 RBM47 RBM47 SAMD4A SAMD4A 3 RBFOX1 RBFOX1 RBM45 RBM45 2 RBM28 RBM28 density 24 h density 2 FMR1 FMR1 8h CPEB2 CPEB2 QKI QKI 1 MBNL1 MBNL1 1 ENOX1 ENOX1 HNRNPLL HNRNPLL PTBP1 PTBP1 0 0 SRSF7 SRSF7 Proteomics RNA seq −1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0 RBM3 RBM3 HNRNPH2 HNRNPH2 ΔPSI ΔPSI PCBP2 PCBP2 SRSF2 SRSF2 KHDRBS1 KHDRBS1 d EndoC-βH1 cells Human Islets TIA1 PABPN1 TIA1 PABPN1 DAZAP1 DAZAP1 Inclusion TARDBP TARDBP ELAVL1 ELAVL1 SRSF6 SRSF6 SRSF9 SRSF9 PCBP1 PCBP1 SNRNP70 SNRNP70 Exclusion HNRNPA1 HNRNPA1 CELF5 CELF5 PABPC1 PABPC1 p = 0.0348 YBX1 YBX1 1.0 1.5 p = 0.0303 IGF2BP3 IGF2BP3 Inclusion/Exclusion Inclusion/Exclusion 0.8 log2 FC log2 FC 1.0 0.6 –1.20 0.00 1.20 –0.50 0.00 0.50 0.4 0.5 i j downstream ExonIntron 0.2 upstream ExonIntron downstream Intron 0.0 0.0 upstream Intron NT IFNα NT IFNα ELAVL1 — TTT[AG][GT]TT Control FMR1 targets e Upregulated SRSF9 — [GT]G[AG][AT]G[GC][AC] 0.00 RBM3 — [AG]A[AGT]AC[GT]A QKI — ACTAAC[ACG] Transcript Protein log2FC p value FDR FMR1 — [GT]GACA[AG]G –0.02 CPEB2 — C[ACT]TTTTT OASL - 001 NP 003724.1 1.135 0.000 0.032 Log2FC HNRNPA1.[AGT]TAGGG[AT] OASL - 002 NP 937856.1 0.724 0.002 0.059 –0.04 ELAVL1 — TTT[AG][GT]TT Downregulated SRSF9 — [GT]G[AG][AT]G[GC][AC] RBM3 — [AG]A[AGT]AC[GT]A –0.06 QKI — ACTAAC[ACG] FMR1 — [GT]GACA[AG]G p = 0.017 CPEB2 — C[ACT]TTTTT –0.08 HNRNPA1.[AGT]TAGGG[AT] –log p value 0.00 2.00 4.00 alternative promoter. The IFNα-induced RMI2−004 upregulation Mining IFNα signatures to identify T1D therapeutic targets. was confirmed using specific primers in both EndoC-βH1 cells Considering the significant overlap observed between gene pro- and human islets (independent samples) (Fig. 6d). These findings files of IFNα-exposed EndoC-βH1 cells and beta cells from T1D support a double mechanism by which IFNα affects human beta individuals (Supplementary Fig. 2d), mining these common sig- cells, i.e. first a massive change in open chromatin regions natures might identify relevant T1D therapeutic targets. First, the followed by later changes in gene expression and AS (see above) top 150 commonly upregulated genes detected by the RRHO and also AFE usage. analysis of both IFNα-exposed EndoC-βH1 cells and beta cells NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 Fig. 5 Interferon-α changes the alternative splicing landscape. a EndoC-βH1 cells were exposed to IFNα for the indicated time points. The significantly upregulated (red) and downregulated (blue) transcripts were identified using Flux Capacitor (n = 5, |log2FC| > 0.58 and FDR < 0.05). b Frequency of individual alternative splicing events regulated by IFNα (n = 5, |ΔPSI| > 0.2, minimum 5 reads, FDR < 0.05). c Frequency distribution of alternative cassette exon (CEx) events altered by IFNα ((n = 5, ΔPSI) > |0.2| and FDR < 0.05). d Confirmation of the increased exon 4 inclusion in the antiviral gene OASL by IFNα (24 h). cDNA was amplified by RT-PCR using primers located in the upstream and downstream exons of the splicing event and the product evaluated using a Bioanalyzer 2100 (n = 4 (EndoC) and n = 7 (human islets), two-sided paired t-test). e The log2FCs of the proteins coding for OASL-001 and −002 isoforms from IFNα-treated EndoC-βH1 cells proteomics (24 h) (n = 4). f Frequency distribution of retained intron (RI) events altered by IFNα (n = 5, | ΔPSI| > 0.2 and FDR < 0.05). g The protein log2FC values obtained by proteomics analysis of EndoC-βH1 cells exposed to IFNα for 24 h were classified in three categories according to the levels of retained intron ΔPSI (n = 5, mean ± SEM, ANOVA with Bonferroni correction). h Expression of RNA-binding proteins (left) that are significantly modified at mRNA level (FDR < 0.05) after exposure to IFNα and their respective proteins (right) in the indicated time points (n = 4–5). i Positional enrichment of motifs from significantly modified RBPs among regions involved in the regulation of modified cassette exons (CEx) after exposure to IFNα for 24 h. (n = 5, |ΔPSI| > 0.2, FDR < 0.05). j Comparison between the log2FC of a curated list (Supplementary Table 1) of known FMR1 target proteins against the log2FC of the remaining proteins detected by the proteomics of EndoC-βH1 cells exposed to IFNα for 24 h (n = 4, mean ± SEM; two-sided unpaired t-test). Source data are provided as a Source Data file. from T1D individuals were selected (Supplementary Fig. 2d and human and rodent beta cell responses to stressful stimuli53,54, all Fig. 7a) to query the Connectivity Map database49. We focused in experiments were performed in clonal or primary human beta opposite signatures of perturbagens that may reverse the effects of cells/islets. This approach identified very rapid and broad beta cell IFNα. To decrease off-target findings based on individual com- responses to IFNα including: (1) major early modifications in pounds, the analysis was performed considering only the classes chromatin remodeling, which activates regulatory elements; (2) of perturbagens. Four main classes, including bromodomain the key TFs regulating signaling, and the crosstalk between beta inhibitors, potentially reversed the signature from our query (tau cells and immune cells; (3) the functional modules of genes and score < −90) (Fig. 7b). Comparable results were obtained when their regulatory networks; and (4) alternative splicing and first analyzing the intersection of IFNα-exposed pancreatic human exon usage as important drivers of transcript diversity. Finally, an islets and beta cells from T1D individuals (Supplementary integrative analysis led to the identification of two compound Fig. 13a). Bromodomain inhibitors have been shown to prevent classes that reverse all or part of these alterations in EndoC-βH1 autoimmune diabetes in animal models50 and the KD of the cells and human islets and may be potential therapeutic targets bromodomain containing 2 gene (BRD2) induced an opposite for future trials in T1D prevention/treatment. signature to our model (Supplementary Fig. 13b). Pre-treatment During viral infection a prompt innate immune response, of EndoC-βH1 cells with two bromodomain inhibitors decreased mediated to a largest extent via type I interferons, is critical to both IFNα-induced HLA-I and CXCL10 induction, with no control virus replication and spreading55. In line with this, changes in CHOP (DDIT3) expression (Fig. 7c, e) or in apoptosis exposure of human beta cells to IFNα leads to changes in chro- induced by IL1β + IFNα (Fig. 7d, f). In human islets, these matin accessibility already at 2 h, which correlates with sub- inhibitors induced a ~30% decrease in IFNα-induced HLA-I sequent changes in mRNA and protein expression at 8 and 24 h. expression and a 90% reduction in CXCL10 expression; at least in The majority of these regions are localized distally to TSSs, the context of I-BET-151, there was a 60% reduction of the ER indicating that they may act primarily as distal regulatory ele- stress marker CHOP (DDIT3) (Supplementary Fig. 13c, d). ments. Interestingly, these regions were enriched in motifs of Next, we searched for clinically approved drugs (DrugBank islets-specific TFs, suggesting that tissue-restricted characteristics 5.151) among the PPI network of the WGCNA module #2 regulate the local responses during insulitis, as we have recently (Fig. 4f), with a view to possible drug repurposing. Module #2 is described for the cytokines IL1β + IFNγ19. This could explain the particularly interesting in this context as it recapitulates many of preferential expression of HLA class I (both the classical ABC the key biological processes induced by IFNα (Fig. 4g), and members and the presently described inhibitory HLA-E) by because ~50% of its members were also present among the most pancreatic islets in comparison to the surrounding exocrine upregulated genes from the RRHO analysis (Supplementary pancreas. Islet HLA class I overexpression is a key finding during Fig. 2d). An interesting target recognized as a hub for different T1D development56, contributing for the recruitment of auto- drugs was the kinase JAK1 (Fig. 8a) and its inhibitor baricitinib, reactive CD8+ T cells that selectively attack beta cells1. IFNα also which has shown promising effects in the treatment of human induces pathways involved in protein modification (ubiquitina- rheumatoid arthritis52. Baricitinib prevented IFNα-induced tion, sumoylation, etc), degradation (proteasome, etc) and ER mRNA expression of HLA-I, CXCL10 and CHOP (DDIT3) in stress, which can generate neoantigens14. EndoC-βH1 cells (Fig. 8b) and human islets (Fig. 8c) and it The IRF and STAT family members are master TFs involved in completely protected EndoC-βH1 cells (Fig. 8d) and human islets IFN-I signaling2. Viruses have developed several species-specific (Fig. 8e) against the pro-apoptotic effects of IFNα + IL1β. mechanisms to antagonize STAT1 and STAT2 activation55. For Furthermore, baricitinib decreased the cell surface protein instances, the NS5 protein of Zika virus degrades human but not expression of MHC class I by >90% in EndoC-βH1 cells (Fig. 9a) mouse STAT257. In the present work, we confirmed the impor- and human islets (Fig. 9b, c). tance of both STAT1 and 2 for INFα signaling in beta cells, and observed that their individual KD is compensated in most cases Discussion by the remaining member, as a possible backup mechanism to We presently modeled the initial changes observed in the islets of protect against pathogens58. Interestingly, IRF1 seems to be a Langerhans during T1D by performing an integrated multi-omics critical regulator of the IFNα-mediated “defense” responses in approach in EndoC-βH1 cells exposed to the early cytokine IFNα. beta cells, including induction of checkpoint proteins such as The model was validated using human islets RNA-seq and PDL1 and HLA-E (present data), and the suppressors of cytokine independent experiments using the same human beta cell line, signaling 1 and 3 (SOCS3) (ref. 59 and present data). This stands pancreatic human islets and FACS-purified human beta cells. Of in contrast to its pro-inflammatory effects in immune cells60. In relevance, taking into account the major differences between line with the possible role for IRF1 in dampening islet 10 NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 ARTICLE a b Chromosome 7 300 300 33.1 mb 33.09 mb 33.08 mb 33.07 mb 33.06 mb NT5C3A-001 NT5C3A-002 Genes 200 200 NT5C3A-004 Events 1500 Ctl 100 100 0 12,000 IFNα 0 0 8h 24 h 8h 24 h 0 10 CGAE TSS 0 c Chromosome 16 EndoC-βH1 cells Human Islets 11.35 mb 11.37 mb 11.39 mb 11.41 mb 11.43 mb RMI2 RMI2 NT5C3A-004 NT5C3A-004 SOCS1 TNP2 PRM1 snoU13 Y_RNA TNP2 MIR548H2 RMI2 SNORA48 PRM3 RMI2 AC009121.1 Beta-actin Beta-actin PRM2 RMI2 PRM2 RMI2 45 p = 0.0002 p = 0.0065 NT5C3A-004/Beta-actin NT5C3A-004/Beta-actin Ctl 0.8 5 0 45 4 0.6 IFNα 3 0 0.4 1500 2 Ctl 0.2 1 0 1500 0.0 0 IFNα NT IFNα NT IFNα 0 CAGE 15 TSS 0 175 POL2 0 d RMI2-004 RMI2 -201 Primer 1 Primer 2 11.406 mb 11.407 mb 11.408 mb 11.409 mb 11.41 mb 11.411 mb RMI2-004 15 Y_RNA RMI2 EndoC-βH1 cells Human Islets Ctl RMI2-004 RMI2-004 0 15 Beta-actin IFNα Beta-actin 0 600 p < 0.0001 p = 0.0007 Ctl 2.0 1.5 RMI2-004/Beta-actin 0 RMI2-004/Beta-actin 600 1.5 IFNα 1.0 0 1.0 CAGE 8 POL2 TSS 0 30 0.5 25 20 0.5 15 10 0 0.0 0.0 NT IFNα NT IFNα Fig. 6 Changes in the alternative transcription start site (TSS) initiation increase the repertoire of IFNα-regulated transcripts. a The tool SEASTAR42 was used to estimate the frequency of differential alternative first exon (AFE) usage induced by IFNα in human beta cells. The total number of IFNα- dependent AFEs events (left) and number of genes with AFEs (right) in the indicated time point are shown (n = 5, ΔPSI > |0.2|, FDR < 0.05). b View of the NT5C3A locus showing the transcripts with AFE usage, the RNA-seq (red) signals of EndoC-βH1 cells exposed or not to IFNα and the CAGE TSSs information (black scale)45 (upper panel). Confirmation of the AFE usage identified by SEASTAR in the gene NT5C3A (lower panel). cDNA was amplified by RT-PCR using primers located in the AFE and in its downstream exon. The PCR products were analyzed by automated electrophoresis using a Bioanalyzer 2100 and quantified by comparison with a loading control. The values were then corrected by the housekeeping gene β-actin. (n = 4 (EndoC) and n = 6 (human islets), two-sided paired t-test). c View of the RMI2 locus showing all the transcripts in this region, the ATAC-seq (blue) and the RNA- seq (red) signals of EndoC-βH1 cells exposed or not to IFNα for 24 h, the CAGE TSSs information (black scale)45 and RNA polymerase II ChIP-seq signal of human K562 cells exposed to IFNα (black)48. A higher magnification of the RMI2-004 locus is presented below (image representative of 4–5 independent experiments). d Confirmation of the AFE usage in the gene RMI2. Genome mapping (upper part) showing the genomic regions used to design-specific primers located in the AFE of the transcript RMI2-004 and in its downstream exon. The PCR product was analyzed by automated electrophoresis using a Bioanalyzer 2100 and quantified by comparison with a loading control. The values were then corrected by the housekeeping gene β-actin. (n = 4 (EndoC) and n = 6 (human islets), two-sided paired t-test). Source data are provided as a Source Data file. inflammation, systemic knockout of IRF1 prevents autoimmune diversity63. Pancreatic beta cells share many characteristics with diabetes in NOD mice61, whereas IRF1 deletion in islets is neuronal cells, including analogous signal transduction, devel- associated with shorter mouse allograft graft function and opmental steps and splicing networks64. Both T1D risk genes65 survival62. and the cytokines IL1β + IFNγ10 modify AS in beta cells. We Alternative splicing (AS) is a species, tissue and context-specific presently identified a preferential alternative exon inclusion after post-transcriptional mechanism that expands the number of IFNα exposure and mapped the potentially involved RBPs, which transcripts originated from the same gene thus increasing protein included the upregulated protein Quaking (QKI). QKI activation NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 a 200 b Down-regulated 150 cell type IFNα treatment –log (p value) 100 summary HCC515 HEPG2 MCF7 50 HA1E HT29 Up-regulated A375 A549 PC3 Up-regulated Down-regulated 0 median tau score Classe of perturbagen Type 1 diabetes –99.32 Leucine rich repeat kinase inhibitor –96.42 Bromodomain Inhibitor –93.64 MDM inhibitor –90.60 Integrin subunits beta LOF Up-regulated Connectivity Map –81.04 Serotonin receptor agonist gene 1 –75.65 GPCR subset GOF gene2 Query tool –74.10 HSP inhibitor –73.37 X linked mental retardation group 1 LOF gene3 ... Class of perturbagens –100.00 100.00 c ETOH I-BET-151 ETOH I-BET-151 ETOH d 1.2 p < 0.0001 1.2 p < 0.0001 I-BET-151 p = 0.0218 1.2 40 ETOH CXCL10/Beta-actin (FOLD) p < 0.0001 p < 0.0001 p < 0.0001 CHOP/Beta-actin (FOLD) HLA-I /Beta-actin (FOLD) p = 0.0171 I-BET-151 Apoptosis (%) 30 0.8 0.8 0.8 20 0.4 0.4 0.4 10 0.0 0.0 0.0 0 NT IFNα NT IFNα NT IFNα NT IFNα + IL1β e JQ– JQ– JQ– f JQ1+ JQ1+ JQ1+ 1.2 p < 0.0001 1.2 p = 0.0309 1.2 40 p = 0.0243 CXCL10/Beta-actin (FOLD) p < 0.0001 p < 0.0001 p = 0.0003 JQ– CHOP/Beta-actin (FOLD) HLA-I/Beta-actin (FOLD) p = 0.0134 JQ1+ Apoptosis (%) 30 0.8 0.8 0.8 20 0.4 0.4 0.4 10 0.0 0.0 0.0 0 NT IFNα NT IFNα NT IFNα NT IFNα + IL1β Fig. 7 Mining the type I interferon signature of pancreatic beta cells for identification of potentially T1D therapeutic targets. a The top 150 upregulated genes identified in Supplementary Fig. 2d were used to query the Connectivity MAP database of cellular signatures49. b Connectivity map classes of perturbagens that promote an opposite signature to the one shared between beta cells of T1D individuals and EndoC-βH1 cells exposed to IFNα (Supplementary Fig. 2d). c, e EndoC-βH1 cells were pretreated for 2 h with the bromodomain inhibitors I-BET-151 (1 μM) (c) or JQ1 + (0.4 μM) (e) and then exposed to IFNα for 24 h. Cells were collected and the mRNA expression for HLA class I (ABC), the chemokine CXCL10 and the ER stress marker CHOP (DDIT3) evaluated. Ethanol (vehicle) and an inactive enantiomer (JQ1−) were used as respective controls for I-BET-151 and JQ1+. (n = 5, mean ± SEM, ANOVA with Bonferroni correction for multiple comparisons). d, f Cell viability after exposure to the combination of cytokines IFNα (2000 U/ml) + IL1β (50 U/ml) in the presence or not of the bromodomain inhibitors (n = 3, mean ± SEM, ANOVA with Bonferroni correction for multiple comparisons). in monocytes promotes extensive changes in AS, favoring their shown protective effects in different animal models of auto- differentiation into pro-inflammatory macrophages66. Further- immunity71, including the diabetes-prone NOD mice50. We have more, QKI binds to the genome of RNA viruses and inhibits their now expanded these findings to human beta cells, showing that replication67. A similar mechanism was recently described for two distinctive BET inhibitors (JQ1+ and I-BET-151) decrease FMR168, another RBP induced by IFNα, which controls protein IFNα-induced responses, including HLA class I and chemokine translation in beta cell (present data). Several other RBPs were overexpression. observed as downregulated by IFNα and identified as potential After binding to its receptor, IFNα promotes phosphorylation regulators of IFNα-induced AS events. Thus, inhibition of of two tyrosine kinases, JAK1 and TYK2, which then trigger the ELAVL1 and HNRNPA1 reproduced IFNα-mediated changes in downstream signaling cascade. Chemical inhibition of JAK1 + exon usage. Different RNA viruses can use both ELAVL169 and JAK2 prevents autoimmune diabetes in NOD mice72 and poly- HNRNPA170 to support their replication, indicating that the morphisms associated with decreased TYK2 function are pro- decreased expression of these proteins may provide an additional tective against human T1D73. We presently observed that IFN-triggered antiviral mechanism. These findings suggest that baricitinib, a JAK1/2 inhibitor recently approved for use in during potentially diabetogenic viral infections, RBPs may have a rheumatoid arthritis by the FDA52, decreased all the three hall- dual role: first as splicing regulators and second as regulators of marks previously identified in islets of T1D individuals and viral replication. initiated by IFNα in human beta cells, namely HLA class I In order to identify novel approaches to protect beta cells in overexpression, ER stress and beta cell apoptosis, supporting its T1D, we analyzed the similarities between beta cell signatures future testing in T1D. from T1D donors and those following IFNα exposure, and In conclusion, we have applied a multi-omics approach to compared the top identified genes/pathways with the Con- study the different levels of gene regulation induced by IFNα in nectivity Map49 and the DrugBank51 database. This identified two EndoC-βH1 cells and pancreatic human islets. This in vitro groups of potential therapeutic agents, namely bromodomain and modeling showed strong correlation with the mRNA profile from JAK inhibitors. Bromodomain (BRD) proteins are components of beta cells of T1D individuals. At the genomic level, early chro- chromatin-remodeling complexes that promote chromatin matin remodeling activated cis-regulatory elements, many of decompaction and transcriptional activation. BET inhibitors have them presenting motifs for islets-specific TFs, providing a 12 NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 ARTICLE a Arseni c t rioxi de Baricit inib JAK1 Barici tinib Tof acitin ib Tofaci tin ib Ru xol itin ib Ruxo lit inib Zinc Ca rfilzom ib Pacli taxel Est ramu stine L-Cy steine Pyrido xal Phosph ate Docetaxel b p < 0.0001 p < 0.0001 p < 0.0001 p < 0.0001 p < 0.0001 p < 0.0001 p < 0.0001 p < 0.0001 CXCL10/Beta-actin (FOLD) p < 0.0001 CHOP/Beta-actin (FOLD) HLA-I/Beta-actin (FOLD) 1.5 1.5 1.5 p = 0.0002 p < 0.0001 p = 0.0037 p < 0.0001 p < 0.0001 p < 0.0001 1.0 1.0 1.0 0.5 0.5 0.5 0.0 0.0 0.0 NT IFNα 100 500 1000 4000 NT IFNα 100 500 1000 4000 NT IFNα 100 500 1000 4000 Baricitinib (μM) Baricitinib (μM) Baricitinib (μM) c DMSO DMSO DMSO Baricitinib Baricitinib Baricitinib 1.2 p = 0.0002 1.2 p < 0.0001 1.2 p = 0.033 CXCL10/Beta-actin (FOLD) p = 0.0003 p < 0.0001 CHOP/Beta-actin (FOLD) HLA-I/Beta-actin (FOLD) 0.8 0.8 0.8 0.4 0.4 0.4 0.0 0.0 0.0 NT IFNα NT IFNα NT IFNα d DMSO e DMSO Baricitinib Baricitinib 30 p < 0.0001 50 p < 0.0001 p < 0.0001 p < 0.0001 40 Apoptosis (%) Apoptosis (%) 20 30 20 10 10 0 0 NT IL1β + IFNα NT IL1β + IFNα Fig. 8 Establishing JAK1 inhibition as protective mechanism against IFNα-mediated inflammation and apoptosis. a The PPI network of module #2 was integrated with the DrugBank repository51 using the CyTargetLinker app78 in Cytoscape. A higher magnification on JAK1 is shown. b EndoC-βH1 cells were pretreated with DMSO (NT) or baricitinib at the indicated concentrations for 2 h. Cells were then left untreated (black bars), or treated with IFNα alone (white bars) without or with the presence of different concentrations of baricitinib (purple scale bars) for 24 h and mRNA expression of HLA class I (ABC), CXCL10 and CHOP (DDIT3) analyzed. The values were normalized by the housekeeping gene β-actin and then by the highest value of each experiment considered as 1 (n = 4, mean ± SEM, ANOVA with Bonferroni correction for multiple comparisons). c Human islets were pretreated with baricitinib (4 μM) or DMSO (vehicle) and then exposed or not to IFNα for 24 h in the presence or not of baricitinib. mRNA expression of HLA class I (ABC), CXCL10 and CHOP (DDIT3) was analyzed and values normalized by the housekeeping gene β-actin and then by the highest value of each experiment considered as 1. (n = 3, mean ± SEM, ANOVA with Bonferroni correction for multiple comparisons). d, e EndoC-βH1 cells (d) and human islets (e) were pretreated with DMSO or baricitinib (4 µM) for 2 h. Subsequently, cells were left untreated or treated with IFNα (2000 U/ml) + IL1β (50 U/ml) in the absence or presence of baricitinib for 24 h. Cell viability was evaluated using nuclear dyes by two independent observers. (d (n = 5), e (n = 4), mean ± SEM, ANOVA with Bonferroni correction for multiple comparisons). Source data are provided as a Source Data file. possible mechanism by which tissue-restricted autoimmune dis- be a source of potential neoantigens. Interestingly, IFNα-exposed eases might arise. Post-translational modifications, alternative human beta cells upregulate co-inhibitory proteins such as PDL1 splicing and first exon usage were induced by IFNα, likely and HLA-E, which may attenuate or delay the autoimmune expanding the repertoire of proteins and transcripts generated by assault. Finally, the present results provide a useful resource for beta cells in response to this inflammatory stimuli. This can also the discovery of compounds that may be used to reverse the NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications 13 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 a b DMSO was performed with two sequential incubation at 4 °C and one wash in FACS buffer DMSO Baricitinib Baricitinib followed each step. Cells were resuspended in FACS buffer, viability dye was added 100 p < 0.0001 1.2 (DAPI) and cells were sorted on a FACSAria III cell sorter (BD Biosciences). HLA-ABC/HOECHST (RFU) p < 0.0001 p = 0.0062 p = 0.0033 Primary (mouse anti-human NTPDase3, hN3-B3S, www.ectonucleotidases-ab. 80 com) and secondary (Alexa Fluor 546 conjugated donkey anti-mouse IgG (A10036, Positive cells (%) MHC Class I 0.8 Thermo-Fisher Scientific)) antibodies were used with the dilutions described in 60 Supplementary Table 5. Data analysis was carried out with FlowJo software 40 (Version 10). 0.4 20 ATAC sequencing processing and analysis. ATAC sequencing was performed in 0 0.0 four independent experiments for each time point (2 and 24 h)75. For ATAC-seq NT IFNα NT IFNα 50,000 EndoC-βH1 cells were exposed or not to IFNα for 2 or 24 h. After that, the cells were harvested, and the nuclei isolated by using 300 μl of cold lysis buffer c HOECHST INSULIN MHC CLASS I MERGE (10 mM Tris–HCl pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% Igepal CA-630). The nuclei pellet was resuspended in a 25 μl transposase reaction mix containing 2 μl of Tn5 transposase per reaction and incubated at 37 °C for 1 h. The tagmented DNA Untreated was isolated using SPRI cleanup beads (Agencourt AMPure XP, Beckman Coulter). For library amplification two sequential 9-cycle PCR were performed (72 °C for 5 min; 98 °C for 30 s; 9 cycles of 98 °C for 10 s, 63 °C for 30 s; and 72 °C for 1 min; and at 4 °C hold). Finally, the DNA library was purified using the MinElute PCR Purification Kit (Qiagen, Venlo, Netherlands). TapeStation and semi-quantitative PCR assays at target positive and negative controls were performed to ensure the quality and estimate the efficiency of the experiment before sequencing. Libraries IFNα 24 h were sequenced single-end on an Illumina HiSeq 2500. Data processing and ana- lysis is described in Supplementary Methods. RNA-sequencing processing and analysis. Total RNA of five independent experiments with EndoC-βH1 cells and six independent preparation of pancreatic Baricitinib 24 h human islets exposed or not to IFNα for different time points was obtained using the RNeasy Mini Kit (Qiagen, Venlo, Netherlands). RNA integrity number (RIN) values were evaluated using the 2100 Bioanalyzer System (Agilent Technologies, Wokingham, UK). All the samples analyzed had RIN values >9. mRNA was obtained from 500 ng of total RNA using oligo (dT)beads, before it was fragmented and randomly primed for reverse transcription followed by second-strand synthesis to generate double-stranded cDNA fragments. The cDNA undergone paired-end Baricitinib + IFNα 24 h repair to convert overhangs into blunt ends. After 3′-monoadenylation and adaptor ligation, cDNAs were purified. Next, cDNA was amplified by PCR using primers specific for the ligated adaptors. (Illumina, Eindhoven, Netherlands). The generated libraries were submitted to quality control before being sequenced on an Illumina HiSeq 2500. RNA-seq data processing and analysis is described in Supplementary Fig. 9 Baricitinib decreases IFNα-mediated MHC class I protein Methods. expression in beta cells. a EndoC-βH1 cells were pretreated with baricitinib (4 μM) or DMSO and then exposed or not to IFNα for 24 h in the presence Proteomics processing and analysis. EndoC-βH1 cells exposed or not to IFNα were extracted using the Metabolite, Protein and Lipid Extraction (MPLEx) or not of baricitinib. MHC class I (ABC) protein expression was measured by approach. A detailed description of the method used for proteomics processing and flow cytometry. The percentage of positive cells was quantified. (n = 6, analysis is provided in Supplementary Methods. mean ± SEM, ANOVA with Bonferroni correction for multiple comparisons). b, c Dispersed human islets were pretreated with baricitinib (4 μM) or Rank–rank hypergeometric overlap (RRHO) analysis. To compare the signature DMSO (vehicle). Next, cells were left untreated, treated with IFNα alone or induced by IFNα with the one present during insulitis in T1D individuals, we with IFNα in the presence of baricitinib for 24 h. MHC class I intensity was performed the RRHO mapping15. For this goal, a full list of log2FC ranked genes quantified in each condition (b) using Fiji software80 and normalized by the from our RNA-seq of EndoC-βH1 cells and human islets (IFNα vs Control, 24 h) were compared against similarly ranked lists of purified primary beta cells obtained HOECHST intensity to correct for the number of cell per area (n = 3, from individuals with T1D16 and T2D17 (T1D/T2D vs non-diabetic). ANOVA with Bonferroni correction for multiple comparisons, RFU (relative In a RRHO map, the hypergeometric p-value for enrichment of k overlapping fluorescence units)). Immunocytochemistry (ICC) analysis (c) of MHC class genes is calculated for all possible threshold pairs for each experiment, generating a I (ABC) (red), insulin (green) and HO (blue) was performed to confirm matrix where the indices are the current rank in each experiment. The log- transformed hypergeometric p-values are then plotted in a heatmap indicating the MHC class I expression in three independent human islet preparations. degree of statistically significant overlap between the two ranked lists in that Scale bar 10 μm. position of the map. Multiple correction was applied using the Benjamini–Yekutieli FDR correction. effects of IFNα on human pancreatic beta cells, paving the way for potential T1D interventional trials. Dynamic regulatory events miner (DREM) modeling. For reconstructing dynamic regulatory networks, we have used the DREM method23, which integrates times series and static data using an Input-Output Hidden Markov Model Methods (IOHMM), where the TF-DNA interaction information obtained from ChIP-seq Culture of EndoC-βH1 cells and human islets, cell treatment. The human experiments48 was used as the input and our RNA-seq time series expression data pancreatic beta cell line EndoC-βH1 was kindly provided by Dr. R. Scharfmann, as the output. A detailed description of DREM-based modeling is provided in University of Paris, France74. Human islet isolation from 20 non-diabetic organ Supplementary Methods. donors (Supplementary Table 2) was performed in accordance with the local Ethical Committee in Pisa, Italy. The use of pancreatic human islets for this project was approved by the Comité d’Ethique hospitalo-facultaire Erasme-ULB. These Weighted gene co-expression network analysis (WGCNA). On each dataset cells were maintained in culture and treated as described in Supplementary (RNA-seq and proteomics), we obtained modules of genes/proteins of similar Methods. expression profiles using WGCNA28. The soft threshold parameter for the RNA- seq dataset was chosen to be 10 (value to approximate a scale-free topology). Similar parameters were used for the analysis of RNA-seq of pancreatic human FACS-purified human beta cells isolation and treatment. Whole pancreatic islets exposed or not to IFNα. Regarding the proteomics dataset, in order to achieve human islets were exposed or not to IFNα for 24 h. After this period, the islets were an approximated scale-free topology, we first normalized each protein expression dispersed into single cells and surface staining was carried out in FACS buffer (PBS in each temporal group (subtraction by mean and division by standard deviation), with 0.5% BSA and EDTA 2 mM final concentration). Indirect antibody labeling and then selected the soft threshold parameter as 14. After merging the modules 14 NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 ARTICLE using a dissimilarity threshold of 0.25, we identified 8 modules in the RNA-seq Real-time PCR analysis. After harvesting of the cells, Poly(A) + mRNA was dataset and 7 modules in the proteomics dataset. obtained using the Dynabeads mRNA DIRECT kit (Invitrogen) and reverse To analyze module quality, we have used a set of statistics (density and transcribed. Detailed description is provided in Supplementary Methods. separability metrics) from the modulePreservation function of the R package WGCNA29. For this purpose, we resampled the dataset 1000-times to create Western blot, immunocytochemistry and flow cytometry. Detailed description reference and test sets from the original data and evaluate module preservation, together with additional information on western blot, immunocytochemistry and represented as the Zsummary for each module across the resulting networks. Zsummary flow cytometry analysis is provided in Supplementary Methods. > 2 indicates moderate preservation and Z > 10 high quality/preservation for each module29. To evaluate WGCNA module preservation in independent samples, we used the same R function, but in this case applying metrics based on module Immunofluorescence. After dewaxing and rehydration, samples were subjected to density and intramodular connectivity to give a composite statistic Zsummary. heat-induced epitope retrieval (HIER) in 10 mM citrate buffer pH 6.0, then probed To evaluate the overlap of RNA-seq and proteomics modules, we considered a in a sequential manner with appropriate antibodies as indicated in Supplementary mRNA to be differentially expressed at 8 or 24 h if its absolute fold-change was Table 4. The relevant antigen–antibody complexes were detected using secondary >1.5 and its FDR < 0.05. Regarding the proteomics dataset, we considered a protein antibodies conjugated with fluorescent dyes (Invitrogen, Paisley, U.K). Cell nuclei to be differentially abundant at 8 or 24 h if the t-test p-value was <0.05. We selected were stained with DAPI. After mounting, images were captured with a Leica only the differentially expressed genes/abundant proteins in the identified AF6000 microscope (Leica, Milton Keynes, UK) and processed using the standard WGCNA modules. We then searched for the overlap between the elements of the LASX Leica software platform (Version 1.9.013747). For quantification studies, RNA-seq and proteomics modules and obtained an overlap p-value randomly selected insulin-containing islets (ICIs) from individuals with or without (hypergeometric probability). We retained overlapping modules with a FDR < 0.05 diabetes were imaged, in addition to insulin-deficient islets (IDIs) from individuals and a minimum of 10 common elements. with diabetes. Thirty ICIs were analyzed from 6 independent individuals (5 islets per individual), 20 IDIs were analyzed from 4 independent individuals (5 islets per individual) and 30 ICIs were analyzed from 6 independent control individuals (5 Protein–protein interaction network analysis. The inBio Map protein–protein islets per individual). The mean fluorescence intensity (MFI) arising from detection interaction (PPI) network database31 was obtained from https://www.intomics. of HLA-E was measured using LASX Leica quantification software. com/inbio/. We first restricted the network to contain only the elements expressed in human beta cells based on our RNA-seq database (mean RPKM > 0.5 in at least one condition). For each WGCNA overlapping module, we identified the proteins Therapeutic targets identification. The top 150 upregulated genes shared among in the PPI network with a significantly high number of protein-to-protein con- the RNA-seq of EndoC-βH1 cells and human islets exposed to IFNα for 24 h and nections to the set of elements in the module (FDR < 0.01, and minimum number the RNA-seq of beta cells16 from T1D individuals were identified by the RRHO of connections equal to 5). We considered only networks obtained for the over- analysis. This list of genes was used to query the Connectivity Map dataset of L1000 lapping modules #1, #2 and #5, as the other overlapping modules returned empty cellular signatures, which has transcriptional responses of human cells to different PPI networks. We then obtained PPI networks for each WGCNA overlapping chemical and genetic perturbations, using the CLUE platform (https://clue.io)49. modules, involving the original set of module elements, plus the respective iden- To identify compounds potentially reverting the effects induced by interferons in tified connecting proteins. Communities of interacting proteins were identified beta cells, we have focused on perturbagens promoting signatures that were using the EAGLE algorithm76 with the following parameters: CliqueSize threshold: opposite (negative tau score) to our query list. Only perturbagens having a median 6 and ComplexSize threshold: 2. Wordclouds of each community were generated tau score < −90 were considered for further evaluation. using information from geneRIFs terms. Additionally, aiming at potential repurposing of drugs under clinical investigation for treatment of other pathologies, we have integrated the DrugBank database v5.151 with the PPI network obtained from WGCNA module #2 using the Gene regulatory network analysis. A network of regulatory interactions was CyTargetLinker v4.0.078 within Cytoscape v3.6 to build a biological network obtained from RegNetwoks30 (www.regnetworkweb.org). As in the PPI network, annotated with drugs. we first restricted the network to contain only the elements we found to be The small molecules and drugs pointed out by these two approaches were then expressed in the RNA-seq dataset. Similarly to the PPI network analysis, for each validated in vitro as described above to verify their impact on IFNα-induced WGCNA overlapping modules, we identified regulators with a significantly high upregulation of cytokines/chemokines, ER stress markers, HLA class I and beta cell number of regulatory connections to the set of elements in the module (FDR < 0.01, apoptosis. and minimum number of connections equal to 4). We then obtained regulatory networks for each WGCNA overlapping modules, involving the original set of Cell viability assessment. The cell viability is described in details in Supple- module elements, plus the respective identified regulators. mentary Methods. To create a non-redundant dataset of motifs from regions of gained open chromatin, we used the compareMotifs.pl script from the package HOMER77 to merge motifs with a similarity score threshold of 0.7. The remaining motifs were Statistical analysis. Data of the confirmatory experiments are expressed as means mapped to the gain open chromatin regions using the annotatePeaks.pl script. ± SEM. A significant difference between experimental conditions was assessed by paired t-test, unpaired t-test, one-way or two-ways ANOVA followed by Bonfer- roni correction for multiple comparisons as indicated using the GraphPad Prism Transcription factor motif analysis. Sequence composition analysis of de novo program version 6.0 (www.graphpad.com). Results with p ≤ 0.05 were considered motifs was performed using findMotifGenome.pl from the package HOMER77 statistically significant. with parameters ‘-size given -bits -mask’. The motifs having a p ≤ 10−12 and observed in >3.5% of the targets were chosen for subsequent analysis. All de novo matches having a similarity score to known TF motifs higher than 0.7 are shown in Reporting summary. Further information on research design is available in the tables (Fig. 4e and Supplementary Fig. 3c), or when no match was present over the Nature Research Reporting Summary linked to this article. this threshold, the first hit was elected and its score is presented. Data availability Alternative splicing changes validation. Alternative splicing changes identified All raw and processed ATAC and RNA-sequencing data that support the findings of this from RNA-seq were validated by RT-PCR using specifically designed primers study have been deposited in NCBI Gene Expression Omnibus (GEO) with the primary (Supplementary Table 3). To confirm cassette exons, the primers were adjacent to accession code GSE133221 (subseries are GSE133218: RNA-seq of EndoC-βH1 cells, the predicted splicing event. This approach allowed us to discriminate between GSE148058: RNA-seq of human islets, GSE133219: ATAC-seq of EndoC-βH1 cells). The variants based on their fragment sizes. For alternative first exon usage (AFE) proteomics datasets have been submitted to Pride under identifier number PXD014244 validation, we have designed primers spanning regions that are unique to the (http://www.ebi.ac.uk/pride/archive/projects/PXD014244). The network of regulatory isoform of interest (Fig. 6g), and then normalized the results by the housekeeping interactions can be obtained from RegNetwoks (http://www.regnetworkweb.org/ gene β-actin. cDNA was amplified using MyTaq Red DNA polymerase (Bioline, download/RegulatoryDirections.zip). The DrugBank database v5.1 can be downloaded London, UK), and PCR products were analyzed using an Agilent 2100 Bioanalyzer from: https://www.drugbank.ca/releases/5-1-0/downloads/all-full-database. The inBio system (Agilent Technologies, Wokingham, U.K.). The molarity of each PCR band Map protein–protein interaction (PPI) network database can be obtained from: https:// corresponding to a specific splice variant was quantified using the 2100 Expert www.intomics.com/inbio/api/data/map_public/2016_09_12/inBio_Map_core_2016_ Software (Agilent Technologies, Diegem, Belgium), and used to calculate the ratio 09_12.zip. The CAGE peaks from FANTOM5 database can be obtained on: http:// inclusion/exclusion (SE) or isoform-X/β-actin (AFE). fantom.gsc.riken.jp/5/datafiles/phase2.5/extra/CAGE_peaks/. The Connectivity Map database can be accessed using the CLUE platform (https://clue.io). The RNA polymerase II (POLR2A) ChIP-seq of human K562 cells can be obtained from the Small-RNA interference. Transfection was performed using Lipofectamine ENCODE project (GSM935474, https://www.encodeproject.org/experiments/ RNAiMAX (Invitrogen) as described in Supplementary Methods. After that, the ENCSR000FAX/). The Exon Ontology database can be accessed from: http://fasterdb. cells were kept in culture for a 48 h recovery period and subsequently exposed or ens-lyon.fr/ExonOntology/. The information about T1D risk genes can be found on not to IFNα as indicated. Supplementary Table 3 describes the sequences of siRNAs immunobase (www.immunobase.org) and GWAS catalog (https://www.ebi.ac.uk/gwas/). used in the present study. NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications 15 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 The source data underlying Figs. 2c–m, 3a–e, g, i, 5d, g, j, 6b, d, 7c–f, 8b–e, 9a, b and patients recently diagnosed with type 1 diabetes. Diabetes 65, 3104–3110 Supplementary Figs. 4e–h, 5c–m, 6b, d, 12a–c, 12e–g, 13c–d are provided as a Source (2016). data file. 28. Langfelder, P. & Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9, 559 (2008). Received: 15 August 2019; Accepted: 23 April 2020; 29. Langfelder, P., Luo, R., Oldham, M. C. & Horvath, S. Is my network module preserved and reproducible? PLoS Comput Biol. 7, e1001057 (2011). 30. Liu, Z. P., Wu, C., Miao, H. & Wu, H. RegNetwork: an integrated database of transcriptional and post-transcriptional regulatory networks in human and mouse. Database (Oxford) 2015, https://doi.org/10.1093/database/bav095 (2015). References 31. Li, T. et al. A scored human protein-protein interaction network to catalyze 1. Eizirik, D. L., Colli, M. L. & Ortis, F. The role of inflammation in insulitis and genomic interpretation. Nat. Methods 14, 61–64 (2017). beta-cell loss in type 1 diabetes. Nat. Rev. Endocrinol. 5, 219–226 (2009). 32. Choi, U. Y., Kang, J. S., Hwang, Y. S. & Kim, Y. J. Oligoadenylate synthase-like 2. Crow, M. K., Olferiev, M. & Kirou, K. A. Type I interferons in autoimmune (OASL) proteins: dual functions and associations with diseases. Exp. Mol. disease. Annu Rev. Pathol. 14, 369–393 (2019). Med. 47, e144 (2015). 3. Eizirik, D. L. & Op de Beeck, A. Coxsackievirus and type 1 diabetes mellitus: 33. Marques, J. et al. The p59 oligoadenylate synthetase-like protein possesses the Wolf’s footprints. Trends Endocrinol. Metab. 29, 137–139 (2018). antiviral activity that requires the C-terminal ubiquitin-like domain. J. Gen. 4. Ferreira, R. C. et al. A type I interferon transcriptional signature precedes Virol. 89, 2767–2772 (2008). autoimmunity in children genetically at risk for type 1 diabetes. Diabetes 63, 34. Jacob, A. G. & Smith, C. W. J. Intron retention as a component of regulated 2538–2550 (2014). gene expression programs. Hum. Genet. 136, 1043–1057 (2017). 5. Foulis, A. K., Farquharson, M. A. & Meager, A. Immunoreactive alpha- 35. Ray, D. et al. A compendium of RNA-binding motifs for decoding gene interferon in insulin-secreting beta cells in type 1 diabetes mellitus. Lancet 2, regulation. Nature 499, 172–177 (2013). 1423–1427 (1987). 36. Mukherjee, N. et al. Integrative regulatory mapping indicates that the RNA- 6. Meyer, S. et al. AIRE-deficient patients harbor unique high-affinity disease- binding protein HuR couples pre-mRNA processing and mRNA stability. ameliorating autoantibodies. Cell 166, 582–595 (2016). Mol. Cell 43, 327–339 (2011). 7. Marroqui, L. et al. Interferon-alpha mediates human beta cell HLA class I 37. Bruun, G. H. et al. Global identification of hnRNP A1 binding sites for SSO- overexpression, endoplasmic reticulum stress and apoptosis, three hallmarks based splicing modulation. BMC Biol. 14, 54 (2016). of early human type 1 diabetes. Diabetologia 60, 656–667 (2017). 38. Meertens, L. et al. FHL1 is a major host factor for chikungunya virus infection. 8. Cooper, N. J. et al. Type 1 diabetes genome-wide association analysis with Nature 574, 259–263 (2019). imputation identifies five new risk regions. Preprint at https://doi.org/10.1101/ 39. Maurin, T. & Bardoni, B. Fragile X mental retardation protein: to be or not to 120022 (2017). be a translational enhancer. Front. Mol. Biosci. 5, 113 (2018). 9. Onengut-Gumuscu, S. et al. Fine mapping of type 1 diabetes susceptibility loci 40. Pasciuto, E. & Bagni, C. SnapShot: FMRP mRNA targets and diseases. Cell and evidence for colocalization of causal variants with lymphoid gene 158, 1446–1446.e1441 (2014). enhancers. Nat. Genet. 47, 381–386 (2015). 41. Alasoo, K. et al. Genetic effects on promoter usage are highly context-specific 10. Eizirik, D. L. et al. The human pancreatic islet transcriptome: expression of and contribute to complex traits. Elife 8, e41673 (2019). candidate genes for type 1 diabetes and the impact of pro-inflammatory 42. Qin, Z., Stoilov, P., Zhang, X. & Xing, Y. SEASTAR: systematic evaluation of cytokines. PLoS Genet. 8, e1002552 (2012). alternative transcription start sites in RNA. Nucleic Acids Res. 46, e45 (2018). 11. Op de Beeck, A. & Eizirik, D. L. Viral infections in type 1 diabetes 43. Gel, B. et al. regioneR: an R/Bioconductor package for the association analysis mellitus–why the beta cells? Nat. Rev. Endocrinol. 12, 263–273 (2016). of genomic regions based on permutation tests. Bioinformatics 32, 289–291 12. Schreiber, G. The molecular basis for differential type I interferon signaling. J. (2016). Biol. Chem. 292, 7285–7294 (2017). 44. Al-Haj, L. & Khabar, K. S. A. The intracellular pyrimidine 5′-nucleotidase 13. Huang, X. et al. Interferon expression in the pancreases of patients with type I NT5C3A is a negative epigenetic factor in interferon and cytokine signaling. diabetes. Diabetes 44, 658–664 (1995). Sci. Signal 11, eaal2434 (2018). 14. Gonzalez-Duque, S. et al. Conventional and neo-antigenic peptides presented 45. Consortium, F. et al. A promoter-level mammalian expression atlas. Nature by beta cells are targeted by circulating naive CD8+ T cells in type 1 diabetic 507, 462–470 (2014). and healthy donors. Cell Metab. 28, 946–960.e946 (2018). 46. Tranchevent, L. C. et al. Identification of protein features encoded by 15. Plaisier, S. B., Taschereau, R., Wong, J. A. & Graeber, T. G. Rank-rank alternative exons using exon ontology. Genome Res 27, 1087–1097 (2017). hypergeometric overlap: identification of statistically significant overlap 47. Inshaw, J. R. J., Cutler, A. J., Burren, O. S., Stefana, M. I. & Todd, J. A. between gene-expression signatures. Nucleic Acids Res. 38, e169 (2010). Approaches and advances in the genetic causes of autoimmune disease and 16. Russell, M. A. et al. HLA class II antigen processing and presentation pathway their implications. Nat. Immunol. 19, 674–684 (2018). components demonstrated by transcriptome and protein analyses of islet beta- 48. Consortium, E. P. et al. Identification and analysis of functional elements in cells from donors with type 1 diabetes. Diabetes 68, 988–1001 (2019). 1% of the human genome by the ENCODE pilot project. Nature 447, 799–816 17. Xin, Y. et al. RNA sequencing of single human islet cells reveals type 2 (2007). diabetes genes. Cell Metab. 24, 608–615 (2016). 49. Subramanian, A. et al. A next generation connectivity map: L1000 platform 18. Cnop, M. et al. Mechanisms of pancreatic beta-cell death in type 1 and type 2 and the first 1,000,000 profiles. Cell 171, 1437–1452.e1417 (2017). diabetes: many differences, few similarities. Diabetes 54, S97–S107 (2005). 50. Fu, W. et al. Epigenetic modulation of type-1 diabetes via a dual effect on 19. Ramos-Rodriguez, M. et al. The impact of proinflammatory cytokines on the pancreatic macrophages and beta cells. Elife 3, e04631 (2014). beta-cell regulatory landscape provides insights into the genetics of type 1 51. Wishart, D. S. et al. DrugBank 5.0: a major update to the DrugBank database diabetes. Nat. Genet. 51, 1588–1595 (2019). for 2018. Nucleic Acids Res. 46, D1074–D1082 (2018). 20. Eizirik, D. L. & Cnop, M. ER stress in pancreatic beta cells: the thin red line 52. Genovese, M. C. et al. Baricitinib in patients with refractory rheumatoid between adaptation and failure. Sci. Signal 3, pe7 (2010). arthritis. N. Engl. J. Med. 374, 1243–1252 (2016). 21. Gurzov, E. N. et al. Pancreatic beta-cells activate a JunB/ATF3-dependent 53. Eizirik, D. L. et al. Major species differences between humans and rodents in survival pathway during inflammation. Oncogene 31, 1723–1732 (2012). the susceptibility to pancreatic beta-cell injury. Proc. Natl Acad. Sci. USA 91, 22. Rondas, D. et al. Citrullinated glucose-regulated protein 78 is an autoantigen 9253–9256 (1994). in type 1 diabetes. Diabetes 64, 573–586 (2015). 54. Brozzi, F. et al. Cytokines induce endoplasmic reticulum stress in human, rat 23. Ernst, J., Vainas, O., Harbison, C. T., Simon, I. & Bar-Joseph, Z. and mouse beta cells via different mechanisms. Diabetologia 58, 2307–2316 Reconstructing dynamic regulatory maps. Mol. Syst. Biol. 3, 74 (2007). (2015). 24. Colli, M. L. et al. PDL1 is expressed in the islets of people with type 1 diabetes 55. Mesev, E. V., LeDesma, R. A. & Ploss, A. Decoding type I and III interferon and is up-regulated by interferons-alpha and-gamma via IRF1 induction. signalling during viral infection. Nat. Microbiol. 4, 914–924 (2019). EBioMedicine 36, 367–375 (2018). 56. Richardson, S. J. et al. Islet cell hyperexpression of HLA class I antigens: a 25. Andre, P. et al. Anti-NKG2A mAb is a checkpoint inhibitor that promotes defining feature in type 1 diabetes. Diabetologia 59, 2448–2458 (2016). anti-tumor immunity by unleashing both T and NK cells. Cell 175, 1731–1743 57. Grant, A. et al. Zika virus targets human STAT2 to inhibit type I interferon e1713 (2018). signaling. Cell Host Microbe 19, 882–890 (2016). 26. Coupel, S. et al. Expression and release of soluble HLA-E is an 58. Blaszczyk, K. et al. The unique role of STAT2 in constitutive and IFN-induced immunoregulatory feature of endothelial cell activation. Blood 109, 2806–2814 transcription and antiviral responses. Cytokine Growth Factor Rev. 29, 71–81 (2007). (2016). 27. Lundberg, M., Krogvold, L., Kuric, E., Dahl-Jorgensen, K. & Skog, O. 59. Moore, F. et al. STAT1 is a master regulator of pancreatic {beta}-cell apoptosis Expression of interferon-stimulated genes in insulitic pancreatic islets of and islet inflammation. J. Biol. Chem. 286, 929–941 (2011). 16 NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16327-0 ARTICLE 60. Langlais, D., Barreiro, L. B. & Gros, P. The macrophage IRF8/IRF1 regulome NIH-NIDDK-HIRN Consortium 1UC4DK104166-01 and Innovate2CureType1 - Dutch is required for protection against infections and is associated with chronic Diabetes Research Fundation (DDRF). D.L.E. and P.M. have received funding from the inflammation. J. Exp. Med. 213, 585–603 (2016). Innovative Medicines Initiative 2 Joint Undertaking under Grant Agreement No. 115797 61. Nakazawa, T. et al. Complete suppression of insulitis and diabetes in NOD (INNODIA). This Joint Undertaking receives support from the Union’s Horizon 2020 mice lacking interferon regulatory factor-1. J. Autoimmun. 17, 119–125 research and innovation programme and “EFPIA”, “JDRF” and “The Leona M. and (2001). Harry B. Helmsley Charitable Trust”. The DiViD study is funded by The South-Eastern 62. Gysemans, C. et al. Interferon regulatory factor-1 is a key transcription factor Norway Regional Health Authority (grant to K.D.-J.), the Novo Nordisk Foundation in murine beta cells under immune attack. Diabetologia 52, 2374–2384 (2009). (grant to K.D.-J.), and through the PEVNET (Persistent Virus Infection in Diabetes 63. Alvelos, M. I., Juan-Mateu, J., Colli, M. L., Turatsinze, J. V. & Eizirik, D. L. Network) Study Group funded by the European Union’s Seventh Framework Pro- When one becomes many-alternative splicing in beta-cell function and failure. gramme (FP7/2007-2013) under grant Agreement Number 261441 PEVNET. Additional Diabetes Obes. Metab. 20, 77–87 (2018). support was from a JDRF Career Development Award (5-CDA-2014-221-A-N) to S.J.R., 64. Juan-Mateu, J. et al. Neuron-enriched RNA-binding proteins regulate a JDRF research grant awarded to the network of Pancreatic Organ Donors – Virus pancreatic beta cell function and survival. J. Biol. Chem. 292, 3466–3480 (nPOD-V) consortium (JDRF 25-2012-516); an MRC Project Grant MR/P010695/1 (2017). awarded to S.J.R. and N.G.M., a studentship grant from the Norman Family Trust 65. Nogueira, T. C. et al. GLIS3, a susceptibility gene for type 1 and type 2 (to S.J.R. and N.G.M.), a Spanish Ministry of Economy and Competiveness (SAF2017- diabetes, modulates pancreatic beta cell apoptosis via regulation of a splice 86242-R to L.P.) and Marató TV3 (201624.10 to L.P.). L.P. is a recipient of a Ramon variant of the BH3-only protein Bim. PLoS Genet. 9, e1003532 (2013). y Cajal contract from the Spanish Ministry of Economy and Competitiveness (RYC- 66. de Bruin, R. G. et al. Quaking promotes monocyte differentiation into pro- 2013-12864) and M.R. is supported by an FI Agència de Gestió d’Ajuts Universitaris i de atherogenic macrophages by controlling pre-mRNA splicing and gene Recerca (AGAUR) PhD fellowship. Part of the work was performed in the Environ- expression. Nat. Commun. 7, 10846 (2016). mental Molecular Sciences Laboratory, a U.S. Department of Energy (DOE) national 67. Liao, K. C. et al. Identification and characterization of host proteins bound to scientific user facility at Pacific Northwest National Laboratory (PNNL) in Richland, dengue virus 3′ UTR reveal an antiviral role for quaking proteins. RNA 24, WA. Battelle operates PNNL for the DOE under contract DE-AC05-76RLO01830. 803–814 (2018). 68. Soto-Acosta, R. et al. Fragile X mental retardation protein is a Zika virus restriction factor that is antagonized by subgenomic flaviviral RNA. Elife 7, Author contributions M.L.C., L.P., T.O.M. and D.L.E. conceived, designed and supervised the experiments. e39023 (2018). M.L.C., M.I.A., J.L.E.H., A.C.B., A.C., H.R.V, E.S.N., N.G.M. and S.J.R. performed and 69. Sokoloski, K. J. et al. Sindbis virus usurps the cellular HuR protein to stabilize analyzed the experiments. L.K., K.D.-J., M.L., M.A.R, J.J.-M., L.M. and P.M. contributed its transcripts and promote productive infections in mammalian and with material and reagents. M.L.C., M.R.R., B.-J.M.W.-R. and J.-V.T performed bioin- mosquito cells. Cell Host Microbe 8, 196–207 (2010). formatic analyses. M.L.C. and D.L.E. wrote the manuscript, and all authors revised it. 70. Lin, J. Y. et al. hnRNP A1 interacts with the 5′ untranslated regions of M.L.C. and D.L.E. are the guarantors of this work and, as such, have full access to all the enterovirus 71 and Sindbis virus RNA and is required for viral replication. J. data in the study and take responsibility for the integrity of the data and the accuracy of Virol. 83, 6106–6114 (2009). the data analysis. 71. Mele, D. A. et al. BET bromodomain inhibition suppresses TH17-mediated pathology. J. Exp. Med. 210, 2181–2190 (2013). 72. Trivedi, P. M. et al. Repurposed JAK1/JAK2 inhibitor reverses established Competing interests autoimmune insulitis in NOD mice. Diabetes 66, 1650–1660 (2017). The authors declare no competing interests. 73. Marroqui, L. et al. TYK2, a candidate gene for type 1 diabetes, modulates apoptosis and the innate immune response in human pancreatic beta-cells. Diabetes 64, 3808–3817 (2015). Additional information 74. Ravassard, P. et al. A genetically engineered human pancreatic beta cell line Supplementary information is available for this paper at https://doi.org/10.1038/s41467- exhibiting glucose-inducible insulin secretion. J. Clin. Invest. 121, 3589–3597 020-16327-0. (2011). 75. Buenrostro, J. D., Giresi, P. G., Zaba, L. C., Chang, H. Y. & Greenleaf, W. J. Correspondence and requests for materials should be addressed to M.L.C. Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nat. Peer review information Nature Communications thanks David Blodgett and the other, Methods 10, 1213–1218 (2013). anonymous, reviewer(s) for their contribution to the peer review of this work. Peer 76. Shen, H., Cheng, X., Cai, K. & Hu, M.-B. Detect overlapping and hierarchical reviewer reports are available. community structure in networks. Phys. A Stat. Mech. Appl. 388, 1706–1712 (2009). Reprints and permission information is available at http://www.nature.com/reprints 77. Heinz, S. et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in identities. Mol. Cell 38, 576–589 (2010). published maps and institutional affiliations. 78. Kutmon, M., Kelder, T., Mandaviya, P., Evelo, C. T. & Coort, S. L. CyTargetLinker: a cytoscape app to integrate regulatory interactions in network analysis. PLoS ONE 8, e82160 (2013). Open Access This article is licensed under a Creative Commons 79. Tan, G. JASPAR2016: data package for JASPAR 2016. R package version Attribution 4.0 International License, which permits use, sharing, 1.12.0, http://jaspar.genereg.net/ (2019). adaptation, distribution and reproduction in any medium or format, as long as you give 80. Schindelin, J. et al. Fiji: an open-source platform for biological-image analysis. appropriate credit to the original author(s) and the source, provide a link to the Creative Nat. Methods 9, 676–682 (2012). Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the Acknowledgements article’s Creative Commons license and your intended use is not permitted by statutory We are grateful to Dr. Andrea A. Schiavo for help in the implementation of the protocol regulation or exceeds the permitted use, you will need to obtain permission directly from for FACS-purified human beta cells isolation and Isabelle Millard, Anyishaï Musuaya, the copyright holder. To view a copy of this license, visit http://creativecommons.org/ Nathalie Pachera and Michaël Pangerl of the ULB Center for Diabetes Research for licenses/by/4.0/. excellent technical support. This work was supported by grants from the Fonds National de la Recherche Scientifique (FNRS), Welbio CR-2015A-06 and CR-2019C-04, Belgium; the Horizon 2020 Program, T2Dsystems (GA667191); the National Institutes of Health, © The Author(s) 2020 NATURE COMMUNICATIONS | (2020)11:2584 | https://doi.org/10.1038/s41467-020-16327-0 | www.nature.com/naturecommunications 17