RNA-Seq Analysis of IL-1B and IL-36 Responses in Epidermal Keratinocytes Identifies a Shared MyD88-Dependent Gene Signature

IL-36 cytokines have recently emerged as mediators of inflammation in autoimmune conditions including psoriasis vulgaris (PsV) and generalized pustular psoriasis (GPP). This study used RNA-seq to profile the transcriptome of primary epidermal keratinocytes (KCs) treated with IL-1B, IL-36A, IL-36B, or IL-36G. We identified some early IL-1B-specific responses (8 h posttreatment), but nearly all late IL-1B responses were replicated by IL-36 cytokines (24 h posttreatment). Type I and II interferon genes exhibited time-dependent response patterns, with early induction (8 h) followed by no response or repression (24 h). Altogether, we identified 225 differentially expressed genes (DEGs) with shared responses to all 4 cytokines at both time points (8 and 24 h). These involved upregulation of ligands (IL1A, IL1B, and IL36G) and activating proteases (CTSS) but also upregulation of inhibitors such as IL1RN and IL36RN. Shared IL-1B/IL-36 DEGs overlapped significantly with genes altered in PsV and GPP skin lesions, as well as genes near GWAS loci linked to autoimmune and autoinflammatory diseases (e.g., PsV, psoriatic arthritis, inflammatory bowel disease, and primary biliary cholangitis). Inactivation of MyD88 adapter protein using CRISPR/Cas9 completely abolished expression responses of such DEGs to IL-1B and IL-36G stimulation. These results provide a global view of IL-1B and IL-36 expression responses in epidermal KCs with fine-scale characterization of time-dependent and cytokine-specific response patterns. Our findings support an important role for IL-1B and IL-36 in autoimmune or autoinflammatory conditions and show that MyD88 adaptor protein mediates shared IL-1B/IL-36 responses.

IL-36 cytokines have recently emerged as mediators of inflammation in autoimmune conditions including psoriasis vulgaris (PsV) and generalized pustular psoriasis (GPP). This study used RNA-seq to profile the transcriptome of primary epidermal keratinocytes (KCs) treated with IL-1B, IL-36A, IL-36B, or IL-36G. We identified some early IL-1Bspecific responses (8 h posttreatment), but nearly all late IL-1B responses were replicated by IL-36 cytokines (24 h posttreatment). Type I and II interferon genes exhibited time-dependent response patterns, with early induction (8 h) followed by no response or repression (24 h). Altogether, we identified 225 differentially expressed genes (DEGs) with shared responses to all 4 cytokines at both time points (8 and 24 h). These involved upregulation of ligands (IL1A, IL1B, and IL36G) and activating proteases (CTSS) but also upregulation of inhibitors such as IL1RN and IL36RN. Shared IL-1B/IL-36 DEGs overlapped significantly with genes altered in PsV and GPP skin lesions, as well as genes near GWAS loci linked to autoimmune and autoinflammatory diseases (e.g., PsV, psoriatic arthritis, inflammatory bowel disease, and primary biliary cholangitis). Inactivation of MyD88 adapter protein using CRISPR/Cas9 completely abolished expression responses of such DEGs to IL-1B and IL-36G stimulation. These results provide a global view of IL-1B and IL-36 expression responses in epidermal KCs with fine-scale characterization of time-dependent and cytokine-specific response patterns. Our findings support an important role for IL-1B and IL-36 in autoimmune or autoinflammatory conditions and show that MyD88 adaptor protein mediates shared IL-1B/IL-36 responses.
Keywords: epidermis, eTs1, il-1, inflammatory bowel disease, keratinocytes, nF-kappaB, rna-seq inTrODUcTiOn The IL-36 cytokines (IL-36A, IL-36B, and IL-36G) are recent additions to the IL-1 family with potent pro-inflammatory effects on epithelial tissues in skin and other organs (1). These cytokines are produced by diverse immune cells, including T cells, dendritic cells, plasma cells, and monocytes, although in skin epidermal keratinocytes (KCs) appear to be the most prominent source (2). The role of IL-36 in skin pathology was first suggested by transgenic mouse models with IL-36A overexpressed in basal KCs, which yielded epidermal thickening and a pro-inflammatory phenotype exacerbated by deletion of the IL-36 receptor antagonist IL-36RN (3). Deletion of IL-36 receptor in mice was additionally protective against full development of imiquimod-induced dermatitis (4), further suggesting that IL-36 cytokines mediate inflammation secondary to disruption of epidermal homeostasis. More recently, evidence has emerged to suggest that IL-36 has an important role in diverse autoimmune conditions, including rheumatoid arthritis, systemic lupus erythematosus, and inflammatory bowel disease (IBD) (5)(6)(7)(8).
Within the intestine, for example, activation of IL-36 signaling appears to facilitate mucosal healing (5), and consistent with this IL36A and IL36G mRNAs are elevated in mucosa of ulcerative colitis patients (6). These findings have intensified research focus on IL-36 ligands and their regulatory proteins, with the hope that anti-IL-36 therapies can provide treatment options for cutaneous and non-cutaneous diseases, similar to the success of other anticytokine biologics demonstrating efficacy in recent years (9).
The emerging evidence implicating IL-36 in cutaneous disease has supported a role in psoriasis vulgaris (PsV) as well as a severe form of PsV known as generalized pustular psoriasis (GPP) (10). Mutations promoting IL-36A activation have been linked to GPP, including IL36RN missense mutations (11,12) and inactivating mutations of the AP1S3 adaptor protein complex 1 (AP-1) subunit (13). By contrast, genome-wide association studies of PsV have thus far not uncovered disease-associated variants near IL-36A, IL-36B, IL-36G, or the antagonist IL-36RN, suggesting that genetic factors predisposing individuals toward PsV do not directly promote IL-36 hyperactivity. This may reflect differences in the immunological basis of PsV compared with GPP, with PsV and GPP mediated primarily by adaptive vs. innate/ autoinflammatory responses, respectively (14). Nonetheless, PsV and GPP can co-occur in the same patient and possess shared pathognomonic features, such as neutrophil infiltration leading to development of pustules (GPP) or Munro's microabscesses (PsV). With regard to the more common PsV, IL-36 expression is abundant in skin lesions and appears to at least provide a psoriatic biomarker, with higher expression in PsV lesions compared with lesions from other inflammatory skin conditions (e.g., eczema) (15). It is notable that in a large scale meta-analysis of skin samples from 237 psoriasis patients, the gene encoding IL-36G (IL36G) was 1 of only 5 genes with elevated expression in lesional skin from all 237 patients, indicating that elevated IL36G expression is a near-universal feature of PsV lesions (16). It was also recently shown that IL-36A is elevated in synovial tissues of patients with psoriatic and rheumatoid arthritis (compared with osteoarthritis), suggesting a possible role for IL-36 in the development of psoriatic arthritis (PsA) as well (17).
In skin, KCs are both a source of IL-36 and important cellular target that can be activated to undergo proliferation and release additional cytokines and chemokines (18). Within the epidermis, key extracellular events preceding IL-36 receptor stimulation have been well studied, although downstream signaling has been mostly studied in transformed cell lines and not primary skin cells (19). IL-36 cytokines are expressed as inactive precursors but are activated by KC-or neutrophil-derived proteases cathepsin S, cathepsin G, elastase, and proteinase-3 (4,20). The IL-36 receptor is a heterodimer that includes IL-1Rrp2 (IL-36R) and IL-1R accessory protein (IL-1RAcP) subunits, with extracellular immunoglobulin and intracellular toll/IL-1 (TIR) domains shared with other receptors from the IL-1 cytokine family (1). The  IL-1RAcP component of the IL-36 receptor is identical to that of  the IL-1 receptor, which may contribute to overlapping effects of  IL-1 and IL-36. Stimulation of the IL-36 receptor is expected to  promote recruitment of the MyD88 adaptor protein to the TIR  domain and activate the JNK, MAPK, and ERK1/2 pathways (19). This triggers a downstream expression response coordinated by key transcription factors, which in various cell types appear to include NF-kappaB and AP-1 (6,19,21,22). It is unclear if these factors are activated by IL-36 in epidermal KCs, although it was shown that IL-36 cytokine expression is positively correlated with NF-kappaB abundance in PsV skin lesions (23).
This study uses RNA-seq to evaluate the gene expression response of primary epidermal KCs to stimulation by cytokines from the IL-1 family (IL-1B, IL-36A, IL-36B, and IL-36G). These ligands exhibit strong amino acid homology and interact with receptors sharing conserved domains and an identical IL-1RAcP subunit (24). The four cytokines were therefore included in this study to characterize potential differences and similarities, and we additionally replicate experiments at two time points (8 and 24 h) to identify early and late responses. We identify IL-1B-and IL-36-responsive genes, characterize known functions of such genes, and investigate relationships to psoriasis (PsV and GPP) based upon skin lesion transcriptome signatures. We further integrate these findings with those from GWA studies to assess whether IL-1B and IL-36 target genes are genetically associated with psoriasis or other autoimmune conditions. Finally, we assess the functional dependence of such genes on MyD88 adaptor protein activity and downstream transcription factors.

rna-seq read Mapping and Quality cTl
Fastq files containing reads were filtered using Cutadapt to remove Illumina adaptor sequences and low quality reads, with a Phred33 quality score cutoff of 30 (-q option) and minimum read length of 20 bp (-minimum-length option) (26). The fastxtoolkit window based quality filter was also applied to remove reads with quality score less than 30 for 50% or more of the read length (27). Pre-and post-filter quality CTL statistics for each sample were calculated with FastQC (28). Reads that did not map to ribosomal sequences were obtained by initially running tophat2 (29) on reads for each sample, with default parameters and GTF file specifying ribosomal read locations only (UCSC hg19/hg38). Unmapped reads from this initial tophat2 run (i.e., reads not aligning to ribosomal RNA sequences) were then used for a second tophat2 run in which reads were mapped using a full GTF file specifying the locations of all human genes (UCSC hg19/hg38). Samtools was used to sort bam files and calculate alignment statistics (30). HTseq was used to tabulate the number of reads mapping to each known human gene, with only reads mapping unambiguously to a given gene included in the tabulation (option -m intersection-strict) (31). Cufflinks was used to calculate normalized FPKM and FPKM confidence intervals for quantification of gene expression (32). Mapping rates and expression profiling efficiency were calculated using the RSeQC and RNA-SeQC software packages (33,34). Raw and processed sequencing data are available from Gene Expression Omnibus under the accession GSE109182.

Dimensionality reduction and global Visualization of cytokine responses
Global visualization of genome-wide cytokine responses was facilitated by several techniques, including hierarchical cluster analysis, principal component (PC) response vectors (35), selforganizing maps (SOMs) (36), and Chernoff faces (37). SOMs were generated using the R Package "kohonen" (function: somgrid), and Chernoff faces were generated via the R package "aplpack" (function: faces). In the current context, Chernoff faces provide an easily interpreted and intuitive multivariate visualization complementing other methods of dimensionality reduction (37). Chernoff face components were scaled in proportion to cytokine FC estimates (cytokine/CTL) for 15  The gene chosen to represent a given cluster was the gene with lowest average Euclidean distance between its response profile and all other genes in the cluster.

identification of Differentially expressed genes (Degs)
Differentially expressed genes were identified from comparisons between cytokine-treated and corresponding CTL cells at each time point (8 and 24 h). Tests for differential expression were performed using protein-coding genes with detectable expression in at least 25% of samples involved in a given comparison. A gene was considered to have detectable expression in a given sample if the count per million mapped reads (cpm) was greater than 0.25 and if the lower limit on the FPKM 95% confidence interval was greater than 0. To evaluate differential expression, raw counts tabulated for each gene were first normalized using the weighted trimmed mean of M-values method (R package: edgeR and function: calcNormFactors) (38). Tests for differential expression were performed by fitting negative binomial generalized linear models (R package: edgeR and function: glmFit) (39), with negative binomial model dispersions estimated using the Cox-Reid-adjusted likelihood method (R package: edgeR and function: estimateGLMTrendedDisp) (40). Negative binomial models were generated using treatment (cytokine-stimulated vs. CTL cells) and cell line (line A, B, and C) as covariates. Differential expression P-values were obtained by comparing the full model (treatment + cell line) and reduced model (cell line only) within the likelihood ratio test framework (R package: edgeR; function: glmLRT) (39). Raw P-values from differential expression analyses were adjusted to CTL the false discovery rate (FDR) using the Benjamini-Hochberg method (41).

Functional analysis of Degs
Differentially expressed genes were analyzed to assess enrichment for KEGG, GO BP, CC, and MF terms using the conditional hypergeometric test implemented in the GOstats R package (function: hyperGTest) (42). Annotated KEGG pathway diagrams were drawn using the pathview R package (43). Enrichment of medical subject heading (MeSH) terms was evaluated using the hypergeometric test implemented in the meshR R package (function: meshHyperGTest) (44). Likewise, enrichment of Disease Ontology categories was evaluated using the DOSE R package (function: enrichDO) (45). Comparisons with ranked gene lists generated from microarray data comparisons involving cytokinetreated KCs or human skin diseases were performed using a database assembled and described in a previous publication (46). The NHGRI-EBI GWAS catalog (release 2017-07-10) was used to MyD88-Dependent IL-36 Responses in KCs Frontiers in Immunology | www.frontiersin.org January 2018 | Volume 9 | Article 80 determine if DEGs were enriched with respect to genes associated with particular human traits or diseases (hypergeometric test) (47). For analysis of genes near loci identified by GWA studies of PsV or PsA, we used the list of loci generated from a recent GWAS meta-analysis of these disease phenotypes (48).
Transcription Factor Binding site analyses DNA sequence regions 5,000 bp upstream of DEGs were evaluated for increased frequency of binding sites associated with human transcription factors or unconventional DNA-binding proteins (uDBPs). This was done by screening a set of 2,935 motifs compiled from multiple sources as described in a previous study (16). Generalized additive logistic regression models were used to evaluate the significance of motif count differences between DEGs and all other KC-expressed genes (49). Human transcription factor superfamily and family annotations were obtained from the TFClass database (50).

Formalin-Fixed Paraffin-embedded (FFPe) Tissue analysis of PsV and gPP skin lesions
Gene expression in PsV and GPP lesions was compared using microarray samples generated from FFPE tissue samples as described in a previous study (10). These data had been generated from the Affymetrix Human Gene 2.1 ST array platform and are available through the Gene Expression Omnibus accession GSE79704. Microarray data files were normalized using robust multichip average (51), and the linear model procedures implemented in the R limma package were used to test for differential expression (52  il-1B and il-36 expression responses are More Disparate at 8 h of Treatment but harmonize after 24 h Global analyses suggested broadly similar expression responses to IL-1B and IL-36, although hierarchical clustering showed that a minority of genes exhibit cytokine-or time-specific responses ( Figure 1A). Responses to IL-36A, IL-36B, and IL-36G varied more following 8 h of treatment but were more consistent after 24 h (Figures 1B,C). Correspondingly, IL-36 expression responses were more similar to those of IL-1B following 24 h of treatment (0.72 ≤ rs ≤ 0.73) compared with 8 h of treatment (0.57 ≤ rs ≤ 0.65) (Figure 1F). Consistent with these trends, SOMs and Chernoff faces reflective of genome-wide response patterns were more similar among cytokines at 24 h of treatment compared with 8 h (Figures 1G,H). These trends demonstrate that expression responses to IL-1B, IL-36A, IL-36B, and IL-36G vary more with early treatment (8 h) but partially harmonize with increased treatment time (24 h).

Type i and ii interferon genes have Time-Dependent il-1B and il-36 responses
We identified between 788 and 1,747 DEGs significantly altered by IL-1B, IL-36A, IL-36B, or IL-36G following 8 or 24 h of treatment (FDR < 0.10 with FC > 2.00 or FC < 0.50; Figure S3  Genes encoding interferon receptor components were upregulated more strongly by IL-1B and IL-36 at 8 h compared with 24 h, e.g., interferon gamma receptor 1 (IFNGR1), interferon gamma receptor 2 (IFNGR2), and interferon alpha and beta receptor subunit 2 (IFNAR2) (Figures 2A-C). Consistent with this, genes with two or more interferon response factor 1 (IRF1) binding sites were more likely to be increased by IL-1B and IL-36 at 8 h of stimulation but not 24 h ( Figure 2D). In addition, type I and II interferon-responsive genes identified by microarray in human KCs were more strongly elevated by IL-1B/IL-36 following 8 h of treatment compared with 24 h (Figures 2E-G). Of 2,500 genes most strongly induced by IFN-g in human KCs, most were more strongly elevated by IL-1B at 8 h compared with 24 h (Figure 2H). These results demonstrate time-dependent IFN responses to IL-1B and IL-36 in KCs, with early IFN gene induction (8 h) that fades or reverses to repression with increased treatment time (24 h).
early il-1B-specific responses are associated with epithelium Development and Mitosis   In panels (a-c), average FPKM (±1 SE) is shown for each group, and asterisks denote significant differences relative to the control (CTL) treatment at the corresponding time point (paired two-sample t-test; n = 2 or 3 per treatment). (D) Genes with interferon response factor 1 (IRF1) binding sites (5 kb upstream region). The middle 50% of FC estimates is shown for genes with 2+ IRF1 binding sites compared with genes with fewer IRF1 sites (magenta font, horizontal axis: P < 0.05, Wilcoxon rank sum test). The IRF1 position weight matrix is shown (right) along with the IRF1 tetrameric structure (bottom right; NCBI structure database). (e) IFN-induced gene signature scores. IFN-induced genes were identified from microarray studies of IFN-treated keratinocytes (left margin), and the average FC for these genes was calculated in IL-1B/IL-36 experiments (bottom margin). Left margin labels indicate the cytokine concentration (in ml), treatment duration, and GEO series accession number. All cytokine experiments were replicated with at least two samples per treatment. (F) Top 30 IFN-g-induced genes (identified from GSE36287). (g) Top 35 INFa-induced genes (identified from GSE36287).
(h) Self-organizing maps (SOMs). The SOM layout was determined only from IFN-g-induced genes (i.e., 2,500 genes most strongly induced by IFN-g, GSE36287). Colors reflect average FC estimates for IFN-g-induced genes assigned to each SOM region (columns 1 and 2 on left). The final column (yellow-blue) displays the mean FC difference for each cytokine with respect to each SOM region (8 h mean FC-24 h mean FC; log2 scale). robust expression responses consistent across Time and shared by all Four cytokine Treatments (il-1B, il-36a, il-36B, and il-36g) IL-1B and IL-36 cytokines interact with different receptors although these share a common IL-1RAcP subunit that may permit cross-stimulation (1, 10). Among IL-36 cytokines, for example, 70-90% of genes regulated by one IL-36 cytokine were correspondingly altered by another at the same time point (Figures 3A-F). Many of these were similarly altered by IL-1B, and it was possible to identify genes with consistent responses to all four cytokines. Following 8 h of treatment, 493 DEGs were elevated, and 133 DEGs were repressed by all 4 cytokines (Figures 3A,D). Similarly, following 24 h of treatment, 297 DEGs were elevated, and 184 DEGs were repressed by all 4 cytokines, respectively (Figures 3A,D). Intersecting these results, we identified 185 DEGs elevated by all 4 cytokines at both time points, along with 40 DEGs repressed by all 4 cytokines at both time points (Figures 3C,F). These 225 DEGs (185 increased + 40 decreased) define a shared response of epidermal KCs to IL-1B, IL-36A, IL-36B, and IL-36G that is consistent at both early and late time points (Figure 4).
Most il-1B and il-36-responsive Degs identified by rna-seq are not identified as Differentially expressed in Previous Microarray analyses The gene expression response to IL-1B was previously evaluated in monolayer KC cultures (55), and the response to IL-36 has been evaluated in reconstituted human epidermis (56). Both studies used microarrays to quantify gene expression (55,56).
In the prior IL-36 study (56), the IL-36 concentration applied to KCs was 5 µg/ml, whereas in this study the concentration was more similar to that observed in human serum (10 ng/ml). Fold changes estimated by RNA-seq in this study showed positive but limited genome-wide correlation to FCs from these prior analyses (rs ≤ 0.26), usually with increased correlations among genes with high expression compared with those with low expression ( Figure  S7A in il-1B and il-36 amplify cytokine/ chemokine Production and genes associated with leukocyte chemotaxis, neutrophil activation, and Mucosal immunity The 185 elevated DEGs included genes encoding cytokines (IL1A, IL1B, IL23A, IL32, and IL36G) and chemokines (CCL20, CXCL1, CXCL2, CXCL5, CXCL6, and CXCL8) ( Figure 4A). Consistent with this, increased DEGs collectively showed significant associations with cytokine activity, chemokine activity, and cytokine receptor interaction ( Figure 5); furthermore, a significant proportion of increased DEGs were localized to the extracellular space suggesting roles as secreted factors (Figure 5E). The 185 DEGs were also associated with gene categories related to leukocyte recruitment   (e.g., granulocyte migration, taxis, leukocyte chemotaxis, neutrophil activation, and chemoattractant activity), mucosal immunity, and epithelial cell signaling (Figure 5). Compared with these findings, the 40 decreased DEGs exhibited weaker and less significant functional associations but were significantly associated with xenobiotic metabolism, positive regulation of apoptosis, development, and transport (Figures 5B,D). Increased and decreased DEGs were both associated with different aspects of bone or joint physiology (increased DEGs: rheumatoid arthritis, osteoclast differentiation, and vitamin D metabolism; decreased DEGs: bone density; Figure 5).
genes induced by il-1B and il-36 are similarly elevated in skin lesions from PsV, gPP, and Other inflammatory skin Diseases Pro-inflammatory effects of IL-1B and IL-36 may contribute to cutaneous inflammation in psoriatic disease (18). Consistent with this possibility, expression responses to IL-1B and IL-36 in epidermal KCs were correlated with those in lesional psoriasis skin (Figures 6A-H). This analysis was performed with comparisons to expression responses in chronic plaque psoriasis (PsV) estimated from a recent RNA-seq meta-analysis (fresh frozen tissues), and with comparison to responses in PsV and GPP from microarray analysis of FFPE tissues (Figures 6A-H) (10). Each lesional skin signature was positively correlated with IL-1B and IL-36 responses in epidermal KCs, but in all cases the correlation was slightly stronger for GPP compared with PsV lesions (FFPE samples) (Figures 6A-H). In agreement with this trend, the 185 DEGs elevated by all 4 cytokines were more strongly enriched among GPP-increased genes compared with PsV-increased genes (FFPE samples; Figures 6I,J,K). Likewise, the 40 DEGs repressed by all 4 cytokines were more strongly enriched among GPP-decreased genes compared with PV-decreased genes (FFPE samples; Figures 6L,M,N). Expression responses to IL-1B and IL-36 were additionally correlated with transcriptional alterations observed in other skin diseases, such as atopic dermatitis (rs = 0.291), eczema (rs = 0.283), and allergic dermatitis (rs = 0.269) ( Figure S8  genes induced by il-1B and il-36 Overlap significantly with genes near gWas loci linked to autoimmune and autoinflammatory conditions (iBD, Psoriasis, and Psa) The NHGRI-EBI GWAS catalog (47) was used to determine if IL-1B/IL-36-increased DEGs overlapped significantly with genes linked to various human diseases or traits (Figure 7A). This identified multiple autoimmune and/or autoinflammatory conditions for which associated genes overlapped significantly with IL-1B/ IL-36-increased DEGs, including IBD/ulcerative colitis, psoriasis/ PsA, systemic sclerosis, primary biliary cholangitis, atopic dermatitis, lupus, and rheumatoid arthritis (Figure 7A). These trends were reinforced by analyses of genes linked to diseases based upon the Disease Ontology database (57), which revealed that increased DEGs were associated with gastrointestinal disease, diseases involving hypersensitivity reactions, and respiratory disease (Figure 7B).
The association between increased DEGs and IBD was stronger than any other observed with respect to the NHGRI-EBI GWAS catalog (Figure 7A), with approximately 12% of IL-1B/IL-36increased DEGs located within 200 kb of an IBD-associated locus ( Figure 7C). Increased DEGs closest to an IBD-associated locus included ETS1, NOD2, and FERMT1 ( Figure 7D). Genes near loci linked to PsV and PsA by GWA studies also overlapped significantly with increased DEGs (Figures 7A,E,G). Increased DEGs closest to a PsV locus included TNFAIP3, ETS1, TNIP, and ZC3H12C (Figure 7F), while increased DEGs closest to a PsA locus included IL23A, TNIP1, and CSF2 ( Figure 7H).

The il-1B and il-36 Transcriptional response rebalances Pathway activation through Positive and negative Feedback
IL-1B and IL-36 treatment of epidermal KCs induced expression of mRNAs encoding the ligands IL1A, IL1B, IL36A, IL36B, and IL36G ( Figure S10A in Supplementary Material), demonstrating that these cytokines have the ability to positively regulate their own expression. We also noted increased expression of cathepsin S (CTSS), which cleaves and activates IL-36 cytokines ( Figure S10A in Supplementary Material), along with elevated expression of proteins encoding receptor subunits (IL1RAP) or downstream signaling elements (IRAK2, MYD88, and TAB2) ( Figure S10A in Supplementary Material). Some transcription responses were associated with inhibitory negative feedback, however, with increased expression of genes encoding receptor antagonists (IL1RN and  IL36RN), anti-inflammatory ligands (IL1F10 and IL37), and a decoy receptor (IL1R2); furthermore, we observed decreased expression of some signaling proteins (TAB1 and IRAK1BP1) ( Figure S10A in Supplementary Material). These expression shifts of IL-1 family-associated genes could in some cases be replicated by other cytokine treatments applied to epidermal KCs (e.g., TNF, IFN-g, and IL-17A; Figures S10B-D in Supplementary Material). Potentially, this was due to induction of IL1B and IL36G mRNAs by these other cytokines, yielding an expression response closely matching that of IL-1B or IL-36 ( Figure S10A in Supplementary Material). The 40 IL-1B/IL-36-decreased DEGs, for example, were almost uniformly repressed by IL-17A treatment of KCs in a prior microarray study ( Figure S10E in Supplementary Material).

nF-KappaB and eTs1 Transcription Factor Binding sites are enriched in regions Upstream of genes induced by il-1B and il-36
A dictionary of 2,934 motifs known to interact with vertebrate transcription factors and uDBPs was screened to identify motifs  The 185 motifs identified with respect to IL-1B/IL-36increased DEGs were frequently associated with 5-GGAA/ TTCC-3 and 5-AGTC/GACT-3 elements recognized by TFs from the immunoglobulin fold or basic domain superfamilies ( Figure S11A in Supplementary Material). There was additionally strong overrepresentation of motifs associated with TFs from the ETS, SOX, GATA, JUN, FOS, and NF-kappaB families ( Figure  S11B in Supplementary Material). Of the 3 most significantly enriched motifs, 2 interacted with NF-kappaB (Figures S11C,D in Supplementary Material) and consistent with this the 185 IL-1B/ IL-36-increased DEGs included 2 genes encoding NF-kappaB subunit proteins (i.e., NFKB2 and RELB; Figure 4A; Figure  S11F in Supplementary Material). Several increased DEGs belonged to the NF-kappaB signaling pathway (KEGG pathway hsa04064), including upstream regulators (IL1B, LYN, and BIR3) and downstream targets regulated transcriptionally by pathway activation (TNFAIP3, NFKBIA, CXCL2, and CXCL8) ( Figure S12 in Supplementary Material). Motifs known to interact with AP-1 were also significantly enriched in regions upstream of increased DEGs (FDR = 1.45e−4) although were not included among the top-ranked motifs ( Figure S11C in Supplementary Material).
The NF-kappaB motif most strongly enriched among increased DEGs featured a 5-GGAA/TTCC-3 half site ( Figure  S11D in Supplementary Material) that was prominent among the complete set of 185 motifs enriched among increased DEG upstream regions ( Figure S11A in Supplementary Material). Notably, this half site resembled an ETS1 binding site included in our dictionary ( Figure S11E in Supplementary Material), which was also significantly enriched among increased DEG upstream regions ( Figure S11C in Supplementary Material). Consistent with this result, ETS1 was among the 185 IL-1B/IL-36-increased DEGs (Figure 3A; Figure S11G in Supplementary Material).
These response patterns were observed only with IL-1B and IL-36G stimulation of WT and KO KCs, with both genotypes showing similar expression responses following stimulation by IFN-g, IL-17A, TNF, and IL-17A + TNF (Figure 8Q).

DiscUssiOn
The cytokines IL-36A, IL-36B, and IL-36G are members of the IL-1 family that have increasingly been linked to autoimmune and autoinflammatory diseases affecting epithelial tissues (2). These diseases include psoriasis (PsV and GPP) but have recently expanded to include acute generalized exanthematous pustulosis and autoimmune blistering diseases (e.g., dermatitis herpetiformis, bullous pemphigoid, and pemphigus vulgaris) (58,59). To better understand the role of IL-36 in these conditions, we used RNA-seq to evaluate effects of IL-1B and IL-36 in primary epidermal KCs. We identified early IL-1B-specific responses involving genes contributing to epidermal development and mitosis, along with time-dependent IL-1B and IL-36 responses of type I and II interferon genes. There was strong overlap among responses to all four cytokines with shared responses of genes near GWAS loci linked to autoimmune/autoinflammatory conditions, includ-  by our study provide points of reference for future analyses, and our findings highlight gene-level connections between autoimmune/autoinflammatory diseases and an MyD88-dependent pathway activated downstream of IL-1B and IL-36.
IL-36 receptor includes an IL-1RAcP subunit that also binds IL-1B and may contribute to overlapping IL-36 and IL-1B responses in epidermal KCs (1). Indeed, most genes regulated by IL-36 cytokines were correspondingly regulated by IL-1B (Figure 3). Following 8 h of cytokine stimulation, it was possible to identify some DEGs regulated only by IL-1B and not IL-36 cytokines ( Figures S5A,B in Supplementary Material), but after 24 h of stimulation almost all IL-1B regulated genes were correspondingly altered by IL-36 cytokines ( Figures S5C,D in Supplementary Material). The highly similar transcriptional responses of the IL-36 cytokines we examined, IL-36A, IL-36B, and IL-36G, are even more remarkable (Figure 3). These findings prompt the question of why IL-1B and three distinct IL-36 cytokines are necessary when their transcriptional fingerprints are so redundant. Although potentially deleterious in the settings of autoimmune and autoinflammatory diseases, IL-1B and IL-36 have critical roles in immune response and pathogen defense (60,61). This was apparent from our data since genes induced by IL-1B and IL-36 were above all else associated with response to bacterium/LPS and granulocyte migration ( Figure 5A). From this perspective, functional redundancy may be advantageous, particularly considering that IL-1B and IL-36 cytokines may not all be expressed within the same tissue. Moreover, whereas IL-1B may be active when secreted, IL-36 cytokines require extracellular processing to gain full activation, which may influence the rate at which downstream responses can develop and contribute to pathogen defense (1). Functional redundancies suggested by our analyses may thus represent an important feature enabling more rapid defense responses, which would be especially critical for epithelial tissues in which IL-36 cytokines are most highly expressed (60,61).
IL-36 cytokines are generally regarded as pro-inflammatory and may act to promote KC and immune cell activation during acute stages of psoriatic lesion development (18). In agreement with this, our findings reveal several ways in which IL-36 amplifies the cytokine network in epidermal KCs, such as upregulation of stimulatory ligands (IL1A, IL1B, IL36A, IL36B, and IL36G), activating proteases (CTSS), and key signaling components (IL1RAP, IRAK2, MYD88, and TAB2). In these respects, the IL-36 signaling pathway exhibits feed-forward dynamics with the capacity to selfamplify following initial stimulation. Upregulation of cathepsin S (CTSS) by IL-1B and all IL-36 cytokines is particularly significant, since cathepsin S was recently identified as an activating protease for IL-36 (62) and it had been observed that cathepsin S is selectively expressed by psoriatic KCs (63). Since other proteases capable of activating IL-36 cytokines are neutrophil-derived (20), CTSS upregulation would provide a mechanism for selfamplification of the IL-36 cascade independent of inflammatory cell infiltration into psoriatic lesions. Despite these effects, IL-1B and IL-36 cannot be regarded as strictly pro-inflammatory, since both cytokines also upregulated anti-inflammatory ligands (IL36RN, IL1RN, IL37, and IL1F10) and downregulated some signaling components (IRAK1BP1 and TAB1). These changes in gene expression may represent compensatory responses that rebalance the IL-36 signaling system to limit excessive or prolonged activation, and their role in psoriatic or other autoimmune diseases may be equally significant, as suggested by the association between GPP and IL-36RN polymorphisms (64). Our findings thus reinforce the idea that IL-36 signaling depends on the balance of several components within an integrated system, with net activation dependent upon the cumulative actions of all components, including extracellular proteases, cytokine ligands, receptor subunits, and downstream signaling elements (24).
The cytokines most convincingly demonstrated to mediate psoriasis plaque formation are those for which inhibition has been effective in patients, such as IL-23, TNF, and IL-17A (9). Although these cytokines lie outside of the IL-1 family, our findings reveal connections to IL-1B and IL-36 signaling in each case, with the IL-1B/IL-36 cascade functioning as either an upstream regulator or downstream target. In the case of IL-23, our results showed that both IL-1B and IL-36 increased expression of IL23A mRNA, demonstrating that anti-IL-23 therapies target a protein positively regulated by the IL-1B/IL-36 signaling cascade. This effect on IL23A expression was additionally notable as a connection between the IL-1B/IL-36 target gene set and those genes linked by GWA studies to PsV and PsA (Figures 7F,H). On the other hand, we identified several cytokine treatments in epidermal KCs leading to a transcriptional response similar to that observed for IL-1B and IL-36 (e.g., IFN-g, TNF, and IL-17A; Figure S10 in Supplementary Material), and further inspection revealed that these treatments upregulate expression of IL1A, IL1B, and IL36G mRNAs ( Figure S10A in Supplementary Material). IL-1B and IL-36 thus function not only as an important regulator of diseasemediating cytokines (i.e., IL23A) but also as a downstream target of such cytokines (i.e., TNF and IL17A). These results are further supported by studies demonstrating co-localization of IL-36A staining with IL-17 and IL-23 in epidermal sections of GPP and acute generalized exanthematous pustulosis (58), correlated expression of IL-36 cytokines, IL-17A, and TNF in psoriasis lesions (65), and an association between IL-36A and IL-17 levels in serum of patients with autoimmune blistering diseases (59).
The shared genetic basis of autoimmune diseases has long been supported by clustering of diverse autoimmune diseases within certain families (i.e., familial autoimmunity) (66). The role of IL-36 in multiple types of autoimmune and autoinflammatory disease is not completely understood, but so far evidence has linked IL-36 to psoriasis (PsV, PsA, and GPP), rheumatoid arthritis, IBD, systemic lupus erythematosus, and Sjögren's syndrome (2,7). In our study, a significant proportion of genes elevated by IL-1B and IL-36 belonged to the KEGG rheumatoid arthritis pathway (e.g., CCL20, CXCL7, and CXCL1; Figure 5). In agreement with this, IL-36A and IL-36B abundance were significantly elevated in serum of rheumatoid arthritis patients compared with healthy CTLs (8), and IL-36A was also elevated in synovial tissues from rheumatoid arthritis patients as compared with osteoarthritis patients (17). When considering genes linked to diseases based only upon GWAS findings, we observed significant overlap between IL-1B/IL-36-increased genes and genes linked to PsV (e.g., TNFAIP3, ETS1, TNIP1, and ZC3H12C) and PsA (e.g., IL23A, TNIP1, and CSF2) (Figure 7). Surprisingly, however, the strongest overall association we observed based upon GWAS findings was between IL-1B/IL-36-increased genes and IBD (Figure 7A), with a total of 23 increased genes linked to IBD based upon GWAS findings, including ETS1, NOD2, FERMT1, and CXCL5 ( Figure 7D). The two main manifestations of IBD include Th1-mediated Crohn's disease and Th2-mediated ulcerative colitis, and of these we observed a much stronger association to GWAS findings for ulcerative colitis (P = 6.7 × 10 −6 ) compared with Crohn's disease (P = 0.16). Consistent with this, recent studies have demonstrated stronger elevation of mRNAs encoding IL-36 cytokines in the inflamed mucosa of ulcerative colitis patients compared with Crohn's disease (6). In addition to these findings, our results demonstrate significant overlap between IL-1B/IL-36 increased DEGs and genes linked by GWA studies to systemic sclerosis and primary biliary cholangitis ( Figure 7A). The role of IL-36 has not been as well studied in these autoimmune conditions, but our findings justify investigations to better understand the significance of IL-36 signaling in these contexts.
IL-36 receptor stimulation is expected to recruit MyD88 adaptor protein and subsequently activate JNK/ERK pathways, with downstream transcriptional responses coordinated by NF-kappaB and AP-1 (6,19,21,22). This model of IL-36 action has largely been established from studies of transformed cell lines rather than primary cells from disease-relevant tissues (19). In primary epidermal KCs, our findings support dependence of IL-36 responses on MyD88 activity, with abrogation of IL-36 responses in the absence of functional MyD88 (Figure 8). Downstream of MyD88, IL-1B, and IL-36 stimulation increased expression of genes encoding NF-kappaB subunits (NFKB2 and RELB), and we observed enrichment of NF-kappaB binding sites upstream of IL-36-increased DEGs. We identified similar enrichment of AP-1 binding sites with respect to IL-36-increased DEGs (data not shown), although genes encoding AP-1 subunits were not included among the 185 genes most robustly upregulated by IL-1B and IL-36 treatment ( Figure 4A). Activation of NF-kappaB has been observed in PsV skin lesions and is positively correlated with abundance of each IL-36 cytokine (23). The role of NF-kappaB in psoriatic disease has been further demonstrated by associations with CARD14 gain-of-function mutations, which lead to NF-kappaB activation and have been linked to PsV, PsA, and GPP (67). It is therefore surprising that siRNA knockdown of NFKB1 (p50) expression did not inhibit IL-36 responses in epidermal KCs, even though expression of IL-1B/IL-36-sensitive genes was altered by NFKB1 inhibition (Figure S14 in Supplementary Material). It is possible that NF-kappaB-dependent responses to IL-36 can proceed even in the setting of decreased NFKB1 expression, in part due to compensatory collateral responses to NFKB1 knockdown (68), such as increased RELA and NFKB2 expression (Figures S14D,K in Supplementary Material) or increased endogenous production of IL-1B and IL-36G to yield stronger pathway stimulation ( Figures S14A,B in Supplementary Material). Alternatively, although we demons trated significant siRNA knockdown of NFKB1 expression ( Figure S13 in Supplementary Material), it is possible that residual protein abundance remained and was sufficient to permit IL-1B/IL-36 responses. Future studies using CRISPR/Cas9 mutagenesis of NFKB1 or other NF-kappaB subunits will therefore be valuable to further dissect out the role of this transcription factor in the IL-1B/IL-36-MyD88 signaling cascade.
The success of anti-cytokine therapies for chronic plaque psoriasis has raised hope that this approach can be extended to other cytokine targets for new drug development (9). Existing biologics have been remarkably effective but still do not resolve symptoms for all patients and for some patients drug effectiveness can decline over time (9). There thus remains a need for continuing research toward new treatments to incorporate into psoriasis patient care. We here demonstrate unique effects of IL-36 in epidermal KCs and functional associations between IL-36 and other cytokines validated as psoriasis drug targets (i.e., TNF, IL-17A, and IL-23). Our findings thus help to support a rationale for further studies to evaluate the translational potential of IL-36 for drug development in the setting of psoriatic disease and with respect to a broad range of autoimmune and/or autoinflammatory conditions (69).

eThics sTaTeMenT
This study was carried out in accordance with the recommendations of the United States National Institutes of Health with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki.