Research ArticleAutoimmunity

Resolving TYK2 locus genotype-to-phenotype differences in autoimmunity

See allHide authors and affiliations

Science Translational Medicine  02 Nov 2016:
Vol. 8, Issue 363, pp. 363ra149
DOI: 10.1126/scitranslmed.aag1974

TYK2’s balancing act

Determining the biological consequences of the thousands of genetic variants that contribute to common diseases for the purpose of improving health care is challenging. Genetic variants that influence autoimmune diseases have been identified in the tyrosine kinase 2 (TYK2) gene, but conflicting evidence regarding their biological impact obscures the therapeutic potential of TYK2. By resolving this conflict, Dendrou et al. have revealed a genetic effect that drives an optimal degree of immune signaling: low enough to be protective against autoimmunity but high enough to prevent immunodeficiency. These findings indicate that TYK2 may be a potential drug target in a number of autoimmune conditions.

Abstract

Thousands of genetic variants have been identified, which contribute to the development of complex diseases, but determining how to elucidate their biological consequences for translation into clinical benefit is challenging. Conflicting evidence regarding the functional impact of genetic variants in the tyrosine kinase 2 (TYK2) gene, which is differentially associated with common autoimmune diseases, currently obscures the potential of TYK2 as a therapeutic target. We aimed to resolve this conflict by performing genetic meta-analysis across disorders; subsequent molecular, cellular, in vivo, and structural functional follow-up; and epidemiological studies. Our data revealed a protective homozygous effect that defined a signaling optimum between autoimmunity and immunodeficiency and identified TYK2 as a potential drug target for certain common autoimmune disorders.

INTRODUCTION

Genome-wide association studies (GWAS) have helped to define the landscape of genetic predisposition to common polygenic diseases, including immune-mediated, metabolic, cardiovascular, and neurological conditions. Elucidation of the biological consequences of disease-associated genetic variation can advance our understanding of the pathophysiological mechanisms that promote these diseases (1). This may allow for the identification of new therapeutic targets. The development of drugs against new targets can be complicated because insufficient efficacy and safety of new compounds have often led to clinical trial suspension (2), but retrospective analysis indicates that drug development programs supported by human genetic evidence are more likely to be successful (3). Therefore, it is conceivable that resolving the functional implications of disease-associated genetic variation can be used prospectively not only to identify potential targets but also to evaluate them through preclinical risk-benefit prediction before drug discovery and trialing (4).

We have investigated the tyrosine kinase 2 (TYK2) gene region based on previous GWAS reports of different, independent association signals within the region across multiple autoimmune disorders (511). We reasoned that such a complex association pattern could provide a framework for cross-comparisons that could be of additional value in dissecting the functional consequences of associated variants. Moreover, TYK2 itself is an attractive immunological candidate for the gene affected by the associated variants because it encodes a nonreceptor tyrosine kinase that is constitutively expressed across immune cell types. TYK2 has been reported to transduce signals downstream of the type I interferon (IFN), gp130, interleukin-10 receptor 2 (IL-10R2), IL-13Rα, and IL-12Rβ1 cytokine receptor families and thus has a pleiotropic role in host responses to infection and in tumor surveillance (12). The index single-nucleotide polymorphisms (SNPs) that have been reported as representative of the association signals in the region include noncoding variants as well as missense SNPs such as rs34536443 and rs12720356, which are located directly in the TYK2 gene (511). Although both these SNPs have been suggested to reduce the function of TYK2 in mediating cytokine signaling (1315), this is inconsistent with rs34536443 conferring protection against the conditions it is associated with, whereas rs12720356 confers protection against some diseases but risk for others. To resolve this conflict and to thereby uncover the predictive power gained from comprehensively studying the TYK2 region, we have used a precision target analysis approach that includes (i) genetic, epigenetic, and transcriptomic meta-analysis; (ii) genotype-dependent in vitro cytokine signaling assessments; (iii) the use of a humanized mouse model; (iv) health record association analysis; and (v) structural investigations.

RESULTS

Joint analysis of genetic associations in the TYK2 region across autoimmune diseases reveals that rs34536443 has a protective effect

Because at least eight different index SNPs in the TYK2 region have been reported to be autoimmune disease–associated (table S1) (511), we sought to evaluate whether these represent completely distinct associations and to determine whether SNPs or more complex haplotypes drive these associations. We enriched our analysis for genetic data from conditions where noncoding variants in the region have been reported as associated and where comprehensive fine-mapping has not been previously performed. The SNP data used were obtained from 8726 ankylosing spondylitis (5), 4017 Crohn’s disease (6), 16,691 multiple sclerosis (8), 2814 psoriasis (9), and 3871 ulcerative colitis (6) patients and 19,738 controls (6, 8), all of European ancestry, who were genotyped on the Immunochip. For fine-mapping, we imputed untyped SNPs using data from the 1000 Genomes Project to infer haplotypes (16, 17), and we then jointly analyzed the patient and control data by multinomial logistic regression (18). The strongest association signal is seen at rs34536443 (P = 3.87 × 10−21), which is highly likely to be the causative variant because this signal is explained by this SNP alone (Fig. 1A and table S2). After controlling for this effect by conditional analysis, a secondary association was observed at rs9797854 (P = 2.51 × 10−12), which is the index SNP representing an associated haplotype that includes a total of 19 SNPs (Fig. 1A and table S2). Controlling for the first two signals revealed a tertiary association at rs12720356 (P = 5.67 × 10−9), which represents a haplotype that includes two other SNPs within the 90% credible set (Fig. 1A and table S2). Further conditioning indicated no additional associations in the TYK2 region, demonstrating that for the diseases analyzed, these three associations explain the previously reported signals represented by numerous different noncoding variants.

Fig. 1. Genetic associations in the TYK2 gene region across autoimmune disorders.

(A) Joint association signal plots generated through multinomial logistic regression analysis. Primary association with the signal at rs34536443 (red). Secondary association with the rs9797854 index SNP (blue). Tertiary association with the rs12720356 index SNP. Scale bars indicate degree of linkage disequilibrium (LD) (r2) relative to each index SNP (green). Triangles indicate SNPs in the associated SNP cluster (90% credible set). The genes in the region are shown below the signal plots. (B) Odds ratios (ORs) for index SNPs for each of the five diseases and 95% confidence intervals (horizontal bars). For rs34536443, a nonadditive effect was observed, and ORs are shown separately for C/G heterozygosity (Het) and C/C homozygosity (Hom). For rs9797854 and rs12720356, the ORs fit an additive model. (C) Summary of associations of the minor alleles at rs34536443, rs9797854, and rs12720356 with psoriasis (Ps), rheumatoid arthritis (RA), systemic lupus erythematosus (SLE), type 1 diabetes (T1D), ankylosing spondylitis (AS), Crohn’s disease (CD), ulcerative colitis (UC), multiple sclerosis (MS), juvenile idiopathic arthritis (JIA), and primary biliary cirrhosis (PBC). The rs34536443 association for ulcerative colitis is shown for the homozygous state; for all other associations, the allelic OR is depicted.

We next assessed the effect on protection against or risk for disease associated with carrying one or two copies of the minor allele at the index SNPs for each of the three association signals. For both rs9797854 and rs12720356, the allelic dosage effect fits an additive model. However, the rs9797854 minor allele only confers protection against multiple sclerosis, because we discovered no evidence of association with the other autoimmune diseases tested. The rs12720356 minor allele was found to confer a previously unreported risk for ankylosing spondylitis, risk for Crohn’s disease and ulcerative colitis, and protection against psoriasis, but was not associated with multiple sclerosis at all (Fig. 1B and table S3A). In contrast, the rs34536443 minor allele is protective across all of the diseases tested, and its dosage effect shows a previously unrecognized significant departure from additivity (P = 5.97 × 10−4) (table S3B), with minor allele homozygosity conferring more than double the protective effect of carrying one such allele (Fig. 1C and table 3C). Notably, for ulcerative colitis, rs34536443-mediated protection is only observable in the homozygous state. Collectively, the findings of our meta-analysis, along with reports of association with other diseases (table S1), pinpoint rs34536443 as the only SNP identified genome-wide to date that protects against 10 different autoimmune conditions (Fig. 1C). This indicates a greater sharing of genetically determined pathophysiological mechanisms across conditions than has been previously recognized (19), and rs34536443 is therefore the prime variant for informing cross-disease therapeutic development.

Distinguishing the functional impact of the TYK2 region association signals suggests effects on different genes

Rs34536443 leads to the substitution of a highly conserved proline residue with an alanine at position 1104 (P1104A) in the kinase domain of the TYK2 protein (fig. S1A) (20). Given the specificity of this molecular effect and noting that rs34536443 was not found to correlate with TYK2 gene expression differences when analyzing expression quantitative trait locus (eQTL) data sets from unstimulated and stimulated immune cell subsets (2124) and tissues (25, 26) (table S4), we considered whether the secondary and tertiary associations have any influence on TYK2 or not. The SNPs in the 90% credible set that constitute the secondary signal haplotype are all located in the SLC44A2 gene (Fig. 1A and table S2), which is a solute carrier family member that has been implicated as an autoantibody target in autoimmune hearing loss (27). These SNPs include missense variants that are unlikely to have a direct influence on TYK2 but also noncoding variants, two of which are potentially regulatory based on their colocalization with epigenetic marks of open chromatin and weak enhancers in immune cell types (table S2). Meta-analysis of eQTL data sets (2126) provided no support for a correlation of the secondary signal with TYK2 transcript expression (table S4), although an effect for cell types or stimulatory conditions not included in the data sets analyzed cannot be ruled out. However, a lack of an impact on TYK2 is in keeping with the disease specificity of this association, which is in sharp contrast to the association pattern of rs34536443.

The tertiary association signal index SNP, rs12720356, results in an I684S substitution in the TYK2 pseudokinase domain (fig. S1A). This domain is thought to regulate the kinase activity of TYK2 (28). As for rs34536443, there is some evidence of a correlation between the I864S substitution and reduced TYK2 activity in overexpression systems and a small number of lymphoblastoid cell lines upon stimulation with type I IFN, although it is unclear whether this might be a direct effect or could be related to reduced expression of the TYK2 protein (13, 15). However, an overall similar functional impact of rs34536443 and rs12720356 is at odds with their cross-disease association pattern. TYK2 impairment as conferred by the minor alleles at each of the two SNPs initially seems plausible based on the protective effect of these alleles against psoriasis, rheumatoid arthritis, systemic lupus erythematosus, and type 1 diabetes, but by extension, this impairment would confer both protection against and risk for ankylosing spondylitis, Crohn’s disease, and ulcerative colitis (Fig. 1C).

To begin to resolve this discrepancy between the genetic associations and the functional correlations, we introduced the rs34536443 and rs12720356 variants into human embryonic kidney (HEK) 293T cells using CRISPR (clustered regularly interspaced palindromic repeats)–Cas9 (CRISPR-associated protein 9)–mediated genomic modification (fig. S1, B and C). This allowed TYK2 expression to be regulated by its endogenous promoter. The modified cells were likely to have a more uniform genomic background compared to lymphoblastoid cell lines, which were derived from different individuals. Upon stimulation with IFN-β, TYK2 phosphorylation was reduced in an rs34536443-dependent fashion, but not by the rs12720356 genotype, and neither SNP had an impact on the expression of the TYK2 protein (fig. S1, D and E). Given that the relative impact of TYK2 impairment depends on the action of the other Janus kinases (JAKs) that it functions in concert with (12), we also assayed the phosphorylation of the TYK2 and JAK substrate, signal transducer and activator of transcription (STAT). Reduced STAT phosphorylation (pSTAT) was observed by the rs34536443, but not the rs12720356, genotype; this reduction was consistent with an additive genetic effect, suggesting that the nonadditive impact of the rs34536443 genotype on disease risk may mainly relate to phenotypic differences downstream of pSTAT. Rs34536443 C/C minor allele homozygous cell lines showed >50% lower pSTAT compared to G/G major allele homozygous lines, although this was not as low as in TYK2 knockout lines (fig. S2). These data do not support an effect of rs12720356 on TYK2 function. However, because the kinase mediates signaling downstream of multiple cytokine receptor families (12), we next compared the effects of rs34536443 and rs12720356 on signaling induced by different cytokines using primary human immune cells.

The primary but not tertiary association signal influences TYK2-mediated cytokine signaling in human immune cells

To resolve which cytokine signaling pathways are specifically influenced by the SNP genotype, we performed stimulations with at least one cytokine for each of the cytokine receptor families that activate TYK2 (12). We hypothesized that type I IFN, IL-10, IL-12, and IL-23 signaling were most likely to be affected because these signaling pathways are impaired in TYK2-immunodeficient patients, who completely lack the TYK2 protein (29), but it was unclear whether a less severe, SNP-dependent impairment of TYK2 function might differentially influence signal transduction downstream of these cytokines. A previous report has suggested that the rs34536443 genotype might also affect IL-6, in addition to type I IFN and IL-10 signaling in T cells. However, this study assessed only a small number of rs34536443 major allele homozygotes and heterozygotes, and no minor allele homozygotes (14). The impact of TYK2 absence or impairment on IL-13 signaling has not been investigated.

For our analyses, human peripheral blood immune cells were obtained from the Oxford BioBank resource because healthy volunteers could be selected and recalled based on their genotype (fig. S3A). Thus, despite rs34536443 minor allele homozygotes having a frequency in the U.K. general population of only ~0.20%, we managed to analyze seven such individuals, representing a 22-fold enrichment for this genotype in our cohort (total n = 162 individuals). Using optimized stimulation conditions that allowed for a maximal response to cytokine (fig. S3, B and C), we assayed the phosphorylation of the different STAT molecules known to mediate signaling downstream of the different cytokines tested. We found ~80% lower pSTAT positivity in the rs34536443 minor compared to the major allele homozygotes in response to type I IFN across all immune cell subsets when assaying STAT1 and the noncanonical type I IFN transducer STAT3 (fig. S4 and Fig. 2, A and B). The influence of the genotype on pSTAT in the primary cells was generally consistent with a nonadditive effect (fig. S4; Fig. 2, A and B; and table S5). The differences in the percentage of pSTAT-positive cells were mirrored by the differences in the total mean fluorescence intensity of pSTAT in stimulated compared to unstimulated cells (tables S5 and S6). No genotype-dependent effects were found after IL-6, IL-10, or IL-13 administration (Fig. 2, C to E; fig. S4A; and table S6), indicating that the other JAKs can compensate for the impairment of TYK2 at the respective receptors for these cytokines.

Fig. 2. Cytokine-induced pSTAT by the rs34536443 genotype in primary human immune cells.

Percentage of pSTAT+ cells by the rs34536443 genotype for CD4+ naïve T cells, CD4+ memory (Mem) T cells, CD8+ naïve T cells, CD8+ memory T cells, CD19+ B cells, and CD14+ monocytes (Mono) upon stimulation with (A) IFN-α, (B) IFN-β, (C) IL-6, and (D) IL-10. CD14+ monocytes were stimulated with (E) IL-13. CD4+ memory T cells and CD8+ memory T cells were stimulated with (F) IL-12 (top graph) and IL-23 (bottom graph). For each cytokine, a schematic of its receptor chains and the TYK2, JAK, and STAT molecules that mediate its signaling are shown to the left of the respective graphs. For (A) to (E), P values provided are for a nonadditive genetic model. For (E) and (F), Wilcoxon matched-pairs rank tests were used. Error bars represent ±SEM.

Although the type I IFN, IL-6, and IL-10 receptors are constitutively expressed across immune cell subsets and the IL-13 receptor is found on monocytes directly ex vivo, detecting pSTAT in response to IL-12 and IL-23 necessitated cellular preactivation, at least in part to facilitate up-regulation of the receptors that bind these cytokines on CD4+ and CD8+ T cells. After preactivation, there was about 70% less IL-12 and IL-23 signaling in CD4+ memory T cells from rs34536443 C/C homozygotes as assayed by pSTAT positivity (Fig. 2F, fig. S4A, and table S6). Consistently, IL-12– and IL-23–dependent differences in CD4+ T helper (TH) type 1 and TH17 signature cytokine production were detected, respectively (fig. S5A). No effect was found for TH2 signature cytokine production (fig. S5A), suggesting that a previous report of an apparent TYK2-related difference in TH2 under nonpolarizing conditions may be indirectly due to TH1 and/or TH17 perturbations (14). Reduced IL-12 and IL-23 signaling was also observed in CD8+ memory T cells (Fig. 2F and table S6). Notably, of the specific subset of genes whose expression was altered in a genotype-specific fashion upon IFN-β stimulation based on RNA sequencing (RNA-Seq) analysis, there was impairment of IL-12Rβ2 transcript up-regulation in CD8+ T cells from the rs34536443 minor compared to the major allele homozygotes (fig. S5B). This indicates the potential for indirect and direct effects of rs34536443 on IL-12 signaling in this cellular compartment. In contrast, for all cytokines and cell types assayed, there were no significant correlations between rs12720356 genotype and pSTAT (tables S7 and S8). The other two SNPs, rs2304257 and rs5030390, that together with rs12720356 constitute the haplotype that drives the tertiary association signal in the TYK2 region are noncoding. Rs2304257 is proximal to the TYK2 promoter region (fig. S6), but because the tertiary association does not correlate with the differences in TYK2 transcript expression and has no impact on pSTAT as a downstream measure of TYK2 expression and activity, we found no evidence of a functional effect of rs2304257 on the TYK2 gene (table S9). However, the tertiary signal does correlate with the expression of the ICAM4, ZGLP1, and ICAM1 genes located in the TYK2 region (21, 22, 24, 26). Because other SNPs that are not autoimmune disease–associated contribute to the variation of ICAM4 and ZGLP1 expression, the correlation between the tertiary signal and the expression of these two genes is irrelevant to disease (table S9). This reveals the genotype-dependent ICAM1 expression difference in CD14+ monocytes—which could influence diapedesis of these cells or immunological synapse formation (30)—as the only identified phenotypic effect solely explained by the tertiary association. This suggests rs5030390 rather than rs12720356 as the likely causative SNP because it localizes to regulatory elements proximal to the ICAM1 promoter (fig. S6) (31, 32). Therefore, of the three genetic signals we assessed in the TYK2 region, rs34536443 alone has a demonstrable impact on TYK2 function, specifically leading to impaired type I IFN, IL-12, and IL-23 signaling.

A TYK2 missense mutation is protective in a humanized mouse model of EAE

One or more of these cytokine signaling pathways have been pathophysiologically implicated through other genetic, cellular, and clinical data in each of the autoimmune conditions that rs34536443 protects against (19, 33, 34), although in multiple sclerosis IFN-β is administered therapeutically (35). To address this apparent inconsistency, we generated humanized mice carrying the targeted, orthologous missense substitution, which in mice is at amino acid position 1124 (P1124A), rather than 1104. Tyk2 Ala1124 homozygous mice had comparable immune profiles and Tyk2 transcript and protein expression relative to Tyk2 Pro1124 homozygotes, although the phosphorylation of Tyk2 was impaired upon stimulation with IFN-β (fig. S7, A to C). As observed with the primary human immune cells, Tyk2 Ala1124 homozygotes also displayed >50% lower pSTAT positivity upon stimulation with type I IFN, IL-12, and IL-23 signaling across cell types (Fig. 3, A to C). These differences in pSTAT can have downstream phenotypic implications: Tyk2 Ala1124 compared to Pro1124 homozygotes showed a 1.7-fold reduction in the capacity of IFN-β to block the IFN-γ–induced up-regulation of major histocompatibility complex class II on peritoneal macrophages (fig. S7D), whereas lower IL-12 and IL-23 signaling leads to reduced TH1 and TH17 signature cytokine production in vitro (fig. S7, E to G).

Fig. 3. Tyk2 Ala1124 homozygous mice have impaired cytokine signaling and are protected against EAE.

pSTAT in splenic immune cells from Tyk2 Pro1124 (red, denoted by P/P; n = 8) and Ala1124 (blue, denoted by A/A; n = 8) homozygous mice after stimulation with (A) IFN-β, (B) IL-12, and (C) IL-23. (D) Experimental autoimmune encephalomyelitis (EAE) disease course shown as mean clinical score after immunization of Pro/Pro (red triangles), Pro/Ala (green squares), and Ala/Ala1124 (blue triangles) mice. *P < 0.05, comparing Pro/Pro1124 and Pro/Ala1124 to Ala/Ala1124 mice. (E) Representative intracellular IL-17A and IFN-γ staining and proportion of (F) IL-17A+, (G) IFN-γ+ IL-17A+, and (H) IFN-γ+ CD4+ memory T cells in draining lymph nodes at peak EAE (or equivalent time point in mice without disease) after immunization with myelin oligodendrocyte glycoprotein (MOG) peptide/complete Freund’s adjuvant (CFA) or phosphate-buffered saline (PBS)/CFA (negative control). (I) Absolute CD4+ T cell numbers in the central nervous system (CNS) at peak EAE. (J) Proportion of MOG-responsive CD4+ T cells in the CNS at peak EAE. (K) Representative intracellular IL-17A and IFN-γ staining and proportion of (L) IL-17A+, (M) IFN-γ+ IL-17A+, and (N) IFN-γ+ CD4+ memory T cells in the CNS at peak EAE after immunization with MOG/CFA or PBS/CFA. P values were calculated using Mann-Whitney U tests. Error bars represent ±SEM.

Upon immunization with MOG peptide in CFA to induce the multiple sclerosis–like disease EAE, heterozygous mice showed decreased disease incidence compared to Pro1124 homozygotes, whereas Ala1124 homozygosity led to a complete protection against EAE (Fig. 3D and fig. S8, A and B). Cytokine knockout (36) and cytokine receptor knockout (37) mouse studies suggest that the observed protection against disease is likely mediated by the diminished IL-23 and, to a lesser extent, IL-12 signaling, particularly through an impact on CD4+ T cells, which are the main pathogenic effectors in the disease model used. Consistent with this, after immunization, peripheral TH17 and TH1 cell impairment was observed on the basis of reduced signature cytokine production by draining lymph node cells from Ala1124 homozygous mice, although total CD4+ T cell numbers were unaffected (Fig. 3, E to H, and fig. S8, C to I). A genotype-dependent difference in cytokine production, assayed directly ex vivo, was also detected in negative control mice immunized with adjuvant alone (Fig. 3, F to H). This indicates that the phenotypic difference is not entirely antigen-dependent—in keeping with the association of rs34536443 across multiple autoimmune diseases in humans—but is instead exaggerated as the immune response to the immunizing antigen develops. In the CNS, comparing Pro1124 homozygous mice at peak disease and Ala1124 homozygotes at an equivalent time point after immunization, the latter had a negligible number of infiltrating CD4+ T cells, and any such cells did not show a pathogenic profile (Fig. 3, I to N, and fig. S9, A to F).

In EAE, abrogated IFN-β signaling does not alter disease incidence but increases severity (38, 39). In the Tyk2 P1124A heterozygous mice that developed disease, severity was not increased relative to Pro1124 homozygotes (fig. S9, G and H), indicating that a stronger Tyk2 impairment is required to mediate a detrimental enough reduction in IFN-β signaling. However, such an effect could not be seen in Ala1124 homozygotes given the lack of EAE development, because the impact of defective IL-23 and IL-12 signaling precedes any beneficial action of IFN-β.

TYK2 P1104A homozygosity defines an optimum between autoimmunity and immunodeficiency

A putative drug-dependent amelioration of human autoimmune disease by targeted TYK2 modulation could be counterbalanced by safety concerns associated with the decreased type I IFN, IL-12, and IL-23 signaling. Extrapolating from patients with genetically driven deficiencies of TYK2 (29) and other components of these cytokine pathways (40, 41), long-term therapeutic inhibition of the kinase could augment the risk of certain severe and recurrent infections. However, the Oxford BioBank rs34536443 minor allele homozygous individuals analyzed were all self-reportedly healthy and showed no immune profile perturbations apart from the impaired cytokine signaling (fig. S10). Moreover, interrogation of medical records from 249 minor allele homozygotes out of a total of 116,732 genotyped individuals of European ancestry from the U.K. Biobank resource revealed no rs34536443-dependent effects on most health measures, including mortality and malignancy, other than the association with autoimmune disease (table S10). Evidence suggestive of a protective effect against self-reported hypothyroidism/myxoedema (P = 3.62 × 10−4), which is commonly of autoimmune etiology, was also found (table S10C). The minor allele homozygotes also showed no increase in their hospitalization due to mycobacterial, bacterial, viral, or fungal infections (with >80% power for detecting an association with such infections and a risk ratio of 6.0 or above) (Fig. 4A). This indicates that rs34536443 minor allele homozygosity may allow for the minimal amount of TYK2-mediated cytokine signaling needed to prevent immunodeficiency.

Fig. 4. Rs34536443 genotype does not lead to immunodeficiency based on the U.K. Biobank health record analysis and cell surface cytokine receptor expression.

(A) Frequency of individuals in the U.K. Biobank hospitalized due to TYK2 immunodeficiency–relevant infections and categorized by rs34536443 genotype. (B) Cell surface type I IFN receptor chain 1 (IFNAR1) expression on HEK293T cell lines genetically modified at the TYK2 locus using CRISPR-Cas9. (C) Cell surface expression of IFNAR1 according to the rs34536443 genotype. Representative surface expression shown for CD4+ naïve T cells (left) and quantification of IFNAR1 expression across subsets (right). (D) Cell surface expression of the IL-12 receptor chains (IL-12Rβ1 and IL-12Rβ2) on preactivated T cells by the rs34536443 genotype. Representative cytokine receptor expression shown for CD4+ memory T cells (left), with corresponding quantification in CD4+ and CD8+ memory T cells (right). (E) Cell surface expression of the IL-23 receptor chains (IL-12Rβ1 and IL-23R) on preactivated T cells according to the rs34536443 genotype. Representative cytokine receptor expression shown for CD8+ memory T cells (left), with corresponding quantification in CD4+ and CD8+ memory T cells (right). The gray histogram (B) and plots (D and E) show isotype control staining. For (B) and (C), P values were calculated using Mann-Whitney U tests. For (D), P values were calculated using Wilcoxon matched-pairs rank tests. Error bars represent ±SEM.

Unlike the rs34536443 minor allele homozygotes, all individuals with complete loss of the TYK2 protein reported to date have an immunodeficient phenotype. This may be not only because there is a complete lack of kinase function in these individuals but also because the absence of the TYK2 protein has a further effect on signaling as cell surface cytokine receptor stability is lowered (29). Regarding cell surface expression of the IFNAR1 chain of the type I IFN receptor, this was observed with our TYK2 knockout cell lines. Conversely, the rs34536443 minor allele homozygous cell lines displayed no decrease in IFNAR1 cell surface expression (Fig. 4B and fig. S11A), which may explain the slightly higher pSTAT1 positivity detected upon stimulation in these cells relative to the knockout cell lines (fig. S2D). The primary immune cells obtained from the rs34536443 minor compared to the major allele homozygous individuals also showed no reduction in cell surface type I IFN, IL-12, or IL-23 receptors, consistent with no genotype-dependent effect on TYK2 protein expression (Fig. 4, C to E, and fig. S11, B to F).

Thus, if the minimal difference in TYK2-mediated cytokine signaling when comparing TYK2 knockout cells to cells with the rs34536443 minor allele homozygous genotype primarily relates to the role of the TYK2 protein in stabilizing cytokine receptors, this implies that the minor allele homozygosity drives a near-complete loss of the signal transduction function of TYK2. To investigate how this effect is conferred at the molecular level, we determined the x-ray crystal structure of the Ala1057 JAK2 kinase domain to 2.66 Å resolution, noting that the TYK2 and JAK kinase domains are structurally highly homologous (42). Superposition of this structure with that of the Pro1057 JAK2 kinase domain (43) demonstrated an overall very similar arrangement of the principal elements of the protein kinase fold, consistent with the more than 80 JAK kinase structures in the Protein Data Bank, but with one exception. Electron density maps indicated an unambiguous and large change starting with the αFG helix and extending to the succeeding αG helix (Fig. 5A and fig. S12). The αFG helix includes the site of mutation and is unique to JAK kinases. The shift is such that the distance between corresponding Cα atoms in these helices is on the order of 15 Å. The generation of this crystal structure required two elements: first, inclusion of the adenosine triphosphate (ATP)–competitive inhibitor TG101209, a molecule that was found, as expected, at the ATP-binding site far from αFG and αG, and second, crystallization from a specific medium unrelated to that previously used for the Pro1057 JAK2/TG101209 structure (43). An additional molecule of TG101209 was unexpectedly found contacting αG and also a neighboring protein molecule provided by the crystal lattice, a finding without precedent among JAK kinase crystal structures. This unintended arrangement prevented straightforward assignment of the structural changes to the P1057A mutation. Comprehensive screening of crystallization conditions did not identify a pair of similar conditions from which both Pro1057 and Ala1057 JAK2 complexes with TG101209 produced crystals, even extending to nucleation with microcrystals in cross-seeding experiments. This distinctly different crystallization behavior of the Pro1057 and Ala1057 JAK2 proteins is consistent with the lack of crystal formation in exhaustive experiments using the Ala1104 TYK2 protein, despite previous success with Pro1104 TYK2 (44). We concluded that the relaxed conformational restraint characteristic of a change from proline to alanine produces a substantial expansion of the low-energy protein conformational space of the kinase domain but that the details of our new structure should be assessed with caution. Nonetheless, it is probable that such an effect on accessible conformations is intrinsic to our observed homozygous Ala1104 TYK2–dependent cytokine signaling impairment, which defines an optimal degree of TYK2 function that protects against autoimmune disease without causing immunodeficiency (Fig. 5B and table S11).

Fig. 5. JAK2 kinase domain Pro→Ala substitution confers conformational change.

(A) Crystal structure of the TYK2-homologous JAK2 Ala1057 kinase domain superimposed onto the JAK2 Pro1057 kinase domain structure [Protein Data Bank (PDB) accession 4JI9]. Both domains were complexed with the JAK2 inhibitor TG101209 (orange molecule) to stabilize the crystal structures. When the alanine residue is present at position 1057, a 15Å shift is observed in the αFG-αG helices (shown in blue) relative to their conformation when the proline residue is present (helices shown in red). (B) Relationship between genetically determined differences in TYK2-relevant cytokine signaling and risk for autoimmune disease versus severe and recurring infections. The kinase domain Pro→Ala substitution may confer a conformational change such that TYK2Ala/Ala mediates an optimal degree of cytokine signaling that may be low enough to prevent autoimmunity but high enough to prevent immunodeficiency. Cyt, cytokine; CytR, cytokine receptor.

DISCUSSION

By resolving conflicting genetic associations and functional correlations in the TYK2 gene region, our study has singled out a highly protective effect against autoimmunity. This arises from minor allele homozygosity at the rs34536443 SNP, which drives a near-complete loss of TYK2 function and consequently impairs type I IFN, IL-12, and IL-23 signaling. The comparison of genetic associations in the region across multiple autoimmune conditions provided the necessary foundation for hypothesis generation regarding the relative role of the associated variants, with implications for drug development. The cross-disease association pattern alone indicated that, contrary to previous suggestions (13, 15), rs34536443 and rs12720356 cannot have fully analogous biological effects despite both being missense variants in the TYK2 gene. Epigenetic and transcriptomic meta-analysis and genotype-specific cytokine signaling assessments confirmed this hypothesis by revealing that compared to the other associations we investigated, rs34536443 is distinguishable as the only SNP with a demonstrable impact on TYK2. Therefore, TYK2 may be a potential target for drug-dependent inhibition for multiple common autoimmune disorders.

Through the elucidation of the relative influence of naturally occurring TYK2 genetic variation on cytokine signaling, we have delineated a genotype-determined gradient of TYK2 activity that may be predictive of the dose-response relationship between drug-dependent targeting of the enzyme and therapeutic outcome. Given the relatively modest protective effect of rs34536443 heterozygosity, a partial drug-induced impairment of TYK2, as might be achieved with broad-spectrum kinase inhibitors, would be therapeutically inadequate. This is especially the case for ulcerative colitis where the minor allele of the SNP confers protection only in its homozygous state. Conversely, on the basis of the consequences of TYK2 deficiency (29), the benefit imparted by a long-term drug-dependent loss of TYK2 or blockade of its binding to cytokine receptors would be offset by certain severe and recurring infections. Existing drugs that block cytokine signaling upstream of TYK2, such as antibodies that target IL-12/IL-23 in psoriasis, are showing therapeutic promise, but they have been linked to serious infections in some instances (45).

The precise mimicking of the impact of rs34536443 minor allele homozygosity could provide maximal efficacy with minimized risk, especially for the treatment of early-stage disease. However, although the combination of genetic and functional follow-up data may help to prospectively predict the therapeutic outcome of autoimmune disease treatment, a limitation of our study is that definitive demonstration of such an optimal effect will ultimately require the generation or identification and subsequent trialing of an appropriately specific pharmaceutical. This may be aided by additional structural analyses to further interrogate the molecular mechanism by which the TYK2 P1104A substitution exerts its effect and how this exactly relates to the conformational change we observed when comparing the Ala1057 and Pro1057 JAK2 kinase domain structures. It should also be noted that the precise degree to which drug-dependent manipulation of TYK2 activity may have some influence on susceptibility to infection remains to be determined. This may vary by geography-dependent variation in infectious disease prevalence and the presence of emerging pathogens, some of which may subvert TYK2-mediated cytokine signaling responses (46).

Initiatives for the generation and collation of large clinically relevant data sets are now under way to accelerate progress toward precision medicine (47), but the use of these data to help improve health care provision to patients with common complex polygenic diseases is a challenge. Our study demonstrates how the integration of cross-disease genetic data with molecular, cellular, structural, and epidemiological information can help to provide preclinical predictions regarding (i) target identification, (ii) indications, (iii) dose responses, and (iv) side effects to thereby facilitate drug development and aid the definition of a therapeutic optimum for these diseases.

MATERIALS AND METHODS

Study design

The study was approved by the National Research Ethics Service Committee South Central, Oxford B (REC reference no. 10/H0605/5;v.2), and fresh peripheral blood samples were obtained from white European volunteers (without self-reported autoimmune or other chronic immunological/inflammatory disorders) from the Oxford BioBank, which is a random, population-based cohort of healthy adults (n > 5000, at the time of analysis). With informed consent, individuals were recalled for blood donation in a balanced, age- and sex-matched fashion based on preselected genotypes of interest. All murine experiments were approved by the University of Oxford Clinical Medicine Ethical Review Committee and licensed under the Animals Scientific Procedures Act of the U.K. Home Office.

Statistical analysis of disease associations

SNP data in the TYK2 locus (chr19, 10,016,000 to 10,928,000; HG19 coordinates) from 8726 ankylosing spondylitis (5), 4017 Crohn’s disease (6), 16,691 multiple sclerosis (8), 2814 psoriasis (9), and 3871 ulcerative colitis (6) patients and 19,738 population controls (6, 8) genotyped on Immunochip were used to localize associations. SNPs with minor allele frequency of <1%, genotype missing rate of >5% in any one batch or >1% overall, and not in Hardy-Weinberg equilibrium in any cohort (P < 10−5) were removed. Samples with excess genotype missing rate (>10% of all SNPs or 2% of clean SNPs), abnormal rate of heterozygosity, duplicated, first- to third-degree relatives of samples already present, and/or non-European ancestry were excluded. A total of 291 SNPs and 55,857 samples passed quality control. For fine-mapping, we inferred haplotypes using SHAPEIT and imputed untyped SNPs with IMPUTE version 2 using as a reference all samples from phase 3 of the 1000 Genomes Project (from IMPUTE2 website, released December 2013); untyped SNPs with an imputation information score of <0.7 were removed, and all imputed SNPs in the 90% credible sets for each association signal had an imputation information score of 0.971 or above. Association analysis was performed with multinomial logistic regression as implemented in the Trinculo tool (18). For imputed SNPs, analysis was performed with dosages to account for genotype imputation uncertainty. For each SNP tested, an additive effect was assumed, and population structure was accounted for by including the first seven principal components as covariates. Conditional analysis for the identification of independent effects was performed stepwise, by including previously identified effects as covariates (rs34536443 for signal 2; rs34536443 and rs9797854 for signal 3), until no independent genome-wide significant association (P < 5 × 10−8) was observed. The 90% credible sets were derived from the likelihood function of the multinomial logistic model (16). For each of the three signals, likelihoods for each SNP under the alternative model were normalized and posterior probabilities were derived from the assumption of one causal SNP per signal. Likelihoods for signal two were estimated conditioning on signal one and likelihoods for signal three were calculated conditioning on signal one and two.

Epigenetic colocalizations and ATAC-Seq

We investigated regulatory consequences of associated SNPs by aligning credible sets with publicly available and in-house epigenetic modification data. We obtained cell subset–specific chromatin segmentation maps for the locus from 111 reference epigenomes from the Roadmap Epigenomics Consortium, and 18 epigenomes from the Blueprint Epigenome Project (www.blueprint-epigenome.eu). In-house chromatin accessibility maps were generated using ATAC-Seq, as previously described (48), to analyze monocytes from three donors. Cells were isolated and stimulated for 0, 2, and 8 hours with lipopolysaccharide (Invivogen) in accordance with Fairfax et al. (23). Paired-end sequencing reads [100 base pairs (bp)] were generated, as previously described (48). Reads were aligned to the University of California, Santa Cruz (UCSC) HG19 reference genome using Bowtie 2 with the parameters -X2000, -I0, --no-discordant, and --no-mixed. Additional processing removed duplicate and nonuniquely aligned reads. Peak-calling was performed using MACS.

eQTL colocalizations

Publicly available eQTL data sets analyzed included the primary data and summary statistics indicated in tables S4 and S9. For each independent association signal, we obtained 95% confidence sets (16). For each analysis, unconditional or conditional, the Bayes factor for a variant is proportional (under weak assumptions) to the posterior probability that it is the causal variant. The 95% confidence set consists of the most significant variants such that 0.95 of the posterior probability is accounted for. For colocalization of susceptibility and eQTL signals, we computed Bayes factors by comparing models where both signals are caused by the same causal variants against models where each signal is caused by distinct variants.

CRISPR-Cas9–modified cell lines

HEK293T cells (American Type Culture Collection) were modified using the CRISPR-Cas9 system as previously described (49). For rs34536443 minor allele introduction, the 5′-CACCGTGGGGGGGCTCTGGCTGG-3′, 5′-CTCTCACCGTGGGGGGGCTCTGG-3′, and 5′-AGCCAGGCCCGCAGCCCCACCGG-3′ guides and the 5′-ccttcggggtgaccctgtatgagctgctgacgcactgtgactccagccagagcgcccccacgGTGAGAGCCAGGCCCGCAGCCCCACCGGGAGTTTGCTAGAGCAATTAGAAAGGCATAGGCTGGGCCAGGC-3′ sense and 5′-GCCTGGCCCAGCCTATGCCTTTCTAATTGCTCTAGCAAACTCCCGGTGGGGCTGCGGGCCTGGCTCTCACcgtgggggcgctctggctggagtcacagtgcgtcagcagctcatacagggtcaccccgaagg-3′ antisense homology-directed repair (HDR) templates were used, whereas the 5′-ATGCGTCAGATGTCTGGTCC-3′ and 5′-GCCTCACAGCAGATTGTTCT-3′ primer pairs were used for validation by Sanger sequencing. For rs12720356, the 5′-GCTGTGTCTTTACAGATATCA-3′ guide was used with the 5′-CAATTCCCCGCCATCCAGAGGGCAGAAGCAGGCAGGTTGCCCCAGAGCAGCTGTGTCTTTACAGatatcatggtgacagagtacgtggagcacggacccctggatgtgtggctgcggagggagcggggccatgtgcccat-3′ sense and 5′-atgggcacatggccccgctccctccgcagccacacatccaggggtccgtgctccacgtactctgtcaccatgatatCTGTAAAGACACAGCTGCTCTGGGGCAACCTGCCTGCTTCTGCCCTCTGGATGGCGGGGAATTG-3′ antisense templates for the major allele and the 5′-CAATTCCCCGCCATCCAGAGGGCAGAAGCAGGCAGGTTGCCCCAGAGCAGCTGTGTCTTTACAGatagcatggtgacagagtacgtggagcacggacccctggatgtgtggctgcggagggagcggggccatgtgcccat-3′ sense and 5′-atgggcacatggccccgctccctccgcagccacacatccaggggtccgtgctccacgtactctgtcaccatgctatCTGTAAAGACACAGCTGCTCTGGGGCAACCTGCCTGCTTCTGCCCTCTGGATGGCGGGGAATTG-3′ antisense HDR templates for the minor allele. The 5′-CTCTGCTGCATTCCGCATTT-3′ and 5′-CACAGGCCACACACCAGGTA-3′ primers were used for validation by Sanger sequencing. For knockout cell lines, the genomic DNA in the region corresponding to the TYK2 N-terminal 4.1 ezrin, radixin and moesin (FERM) domain was targeted because this is the domain with the most naturally occurring TYK2 knockout mutations (29). For these lines, the 5′-ACTGTGGGAGCTGTCGACCG-3′, 5′-GTGGGAGCTGTCGACCGAGG-3′, 5′-GCATGAGTTTGTGAATGACG-3′, and 5′-GAATGACGTGGCATCACTGT-3′ guides were used. Polymerase chain reaction (PCR) and Sanger sequencing confirmation were performed using the 5′-AGATCCCCAGAGATGCAAGC-3′ and 5′-GGAACCTGCGGAAGACGTTC-3′ primers. Before the cells (1 × 106/ml) were stimulated with IFN-β (1000 ng/ml; PeproTech) for 15 min at 37°C, they were serum-starved for 6 hours.

Immunoblotting

Cells were lysed in 50 mM tris-HCl (pH 8.0), 10% (v/v) glycerol, 25 mM EDTA, 150 mM NaCl, 2 mM dithiothreitol, 0.5% NP-40 (IGEPAL CA-630), 25 mM sodium fluoride, 1 mM sodium orthovanadate supplemented with protease inhibitors (Roche), and 0.1% phosphatase inhibitor cocktail 2 (Sigma-Aldrich). Anti-actin (clone C4, Millipore) and anti–N-terminal TYK2 (clone 51, BD Biosciences) with anti-phosphotyrosine (Y1054/Y1055) TYK2 (Cell Signaling Technology) or anti–N-terminal TYK2 (H-135, Santa Cruz Biotechnology) with anti–phosphotyrosine (Y1054/Y1055) TYK2 (Santa Cruz Biotechnology) were used for human and murine protein detection, respectively. Membranes were stained with IRDye 680LT- and 800CW-conjugated secondary antibodies and visualized and quantified using an Odyssey Infrared Imaging System (LI-COR Biosciences).

Genotyping

Rs34536443, rs9797854, and rs12720356 were genotyped using TaqMan SNP Genotyping Assays (Applied Biosystems), and mouse genotyping was performed using a custom-made TaqMan SNP Genotyping Assay (Applied Biosystems) that enabled discrimination of alleles encoding Tyk2 Pro1124 and Ala1124.

Flow and mass cytometry

For cell surface staining, cells were stained for 30 min and fixed with 1× BD CellFIX (BD Biosciences). Phosflow staining was performed according to the BD Biosciences protocols using BD Cytofix Fixation Buffer and BD Phosflow Perm Buffer III. Intracellular cytokine staining was performed using the BD Cytofix/Cytoperm Fixation/Permeabilization kit.

Anti-human antibodies used were as follows: anti-IFNAR1 (Abgent) along with Alexa Fluor 488 goat anti-rabbit immunoglobulin G (IgG) (Invitrogen); allophycocyanin (APC) anti-CD25 (clones M-A251 and 2A3) and fluorescein isothiocyanate (FITC) anti-Lin1 (BD Biosciences); Alexa Fluor 700/APC/phycoerythrin (PE) anti-CD4 (clone RPA-T4), PE-Cy7 anti-CD11b (clone ICRF44), peridinin-chlorophyll-protein complex (PerCP)–Cy5.5 anti-CD14 (clone HCD14), Pacific Blue anti-CD16 (clone 3G8), Alexa Fluor 700 anti-CD27 (clone M-T271), APC anti-CD28 (clone CD28.2), Pacific Blue anti-CD45RA (clone HI100), APC anti-CD56 clone (HCD56), PerCP-Cy5.5 anti-CD57 (clone HNK-1), APC anti-CD123 (clone 6H6), APC-Cy7 anti-CD193 (5E8), APC anti–IFN-γ (clone B27), APC anti-IgD (clone IA6-2), and FITC/PE anti–IL-4 (clone MP4-25D2) (BioLegend); Alexa Fluor 488 anti-CD3 (clone OKT3), PE-Cy7 anti-CD4 (clone RPA-T4), Alexa Fluor 700/FITC/PE anti-CD8 (clone OKT8), Alexa Fluor 700 anti-CD11c (clone 3.9), FITC anti-CD15 (clone HI98), Alexa Fluor 488/Alexa Fluor 700/APC anti-CD19 (clone HIB19), eFluor 450 anti-CD38 (clone HB7), PE-Cy7 anti-CD127 (clone eBioRDR5), PE anti-CD197 (clone 3D12), and FITC/PE anti–IL-17 (clone N49-653) (eBioscience); FITC anti–IL-12Rβ1 (clone 69310), APC anti–IL-12Rβ2 (305719), and PE anti–IL-23R (clone 218213) (R&D Systems).

Anti-mouse antibodies used were as follows: PE anti–I-A/I-E (clone M5/114.15.2) and PE anti–IL-4 (clone 11B11) (BD Biosciences); APC–eFluor 780/FITC anti-CD11b (clone M1/70), APC anti-CD11c (clone N418), Alexa Fluor 700/eFluor 450 anti-CD19 (clone eBio1D3), Alexa Fluor 700 anti-CD3e (clone 145-2c11), Alexa Fluor 700/PE anti-CD4 (clone RM4-5), APC/eFluor 450/PE-Cy7 anti-CD8a (clone 53-6.7), APC anti-CD25 (PC61.5), eFluor 450 anti-CD44 (IM7), FITC anti-CD49d (clone R1-2), FITC/PerCP-Cy5.5 anti-CD62L (clone MEL-14), PE anti-CD122 (clone TM-b1), eFluor 660 anti-Eomes (clone Dan11mag), APC anti–IFN-γ (clone XMG1.2), FITC anti–IL-17A (clone eBio17b7), eFluor 450 anti–Ly-6C (clone HK1.4), PE-Cy7 anti-Ly6G (clone RB6-8C5), and APC/PerCP-Cy5.5 anti-NK1.1 (clone PK136) (eBioscience). For phosflow staining, the following antibodies were used: Alexa Fluor 647 anti-pSTAT1 (clone 4a), Alexa Fluor 488 anti-pSTAT3 (clone 4/P-STAT3), Alexa Fluor 647 anti-pSTAT4 (clone 38/p-Stat4), and Alexa Fluor 647 anti-pSTAT6 (clone 18/P-Stat6) (BD Biosciences).

Isotype control antibodies used were as follows: APC/FITC/PE mouse IgG1, PE mouse IgG2b (R&D Systems), PE rat IgG2a (BD Biosciences), and rabbit polyclonal IgG (Abcam). Mass cytometry was performed using peripheral blood mononuclear cells (PBMCs) and the Maxpar Human Peripheral Blood Phenotyping Panel Kit (Fluidigm). Data obtained were analyzed using FlowJo (Tree Star). For fluorescence-activated cell sorting (FACS), murine CNS mononuclear cells were first stained for 40 min at room temperature with APC-labeled 1-Ab/MOG38–49 tetramer (4 μg/ml) or negative control 1-Ab/hCLIP tetramer (National Institutes of Health Tetramer Core Facility, Emory University). For the final 20 min, cells were stained at 4°C with the tetramers, Fc block (BD Biosciences) and PE anti-CD4 (clone RM4-5), Alexa Fluor 488 anti-CD11b (clone M1/70), and eFluor 450 anti-CD45.2 (clone 104) (eBioscience).

Human blood cell subset isolation and culturing

PBMCs were separated on Lymphoprep (Axis-Shield). Specific subsets (CD3+ T cells, CD4+ T cells, CD8+ T cells, CD19+ B cells, and CD14+ monocytes) were isolated using respective positive selection MicroBeads (Miltenyi), and for RNA-Seq, cells were serum-starved for 4 hours and stimulated for 12 hours with IFN-β (400 ng/ml). For cytokine stimulations, IFN-α, IFN-β, IL-6, IL-10, or IL-13 (400 ng/ml; PeproTech) was used to stimulate PBMCs (1 × 106 cells/ml) for 15 min in RPMI at 37°C. For IL-12 stimulation, purified CD3+ T cells were serum-starved for 4 hours, activated for 72 hours with phytohemagglutinin (10 μg/ml; Sigma-Aldrich) in RPMI supplemented with 10% human AB serum, 2 mM glutamine, and penicillin/streptomycin (all from Sigma-Aldrich), then stimulated with IL-2 (100 U/ml; PeproTech) for 24 hours, and serum-starved for 4 hours before stimulation with IL-12 (400 ng/ml; PeproTech) for 15 min. For IL-23 stimulation, PBMCs were activated with Dynabeads Human T-Activator CD3/CD28 (Invitrogen) according to the instructions but with a 1:2 bead-to-cell ratio for 5 days, serum-starved for 4 hours, and then stimulated with IL-23 (400 ng/ml; R&D Systems) for 15 min.

RNA sequencing

Human polyA sequencing libraries were prepared according to the Illumina TruSeq protocol. Libraries were sequenced using Illumina HiSeq 4000 (75-bp paired-end reads). Reads were mapped to UCSC HG19 human using TopHat. Uniquely aligned reads were extracted, and read duplicates were removed using Picard (http://picard.sourceforge.net/). HT-Seq was used to quantify gene expression using the gene annotations provided by Illumina’s iGenomes Project. Differential expression was assessed with the R Bioconductor package DESeq2. Gene expression counts were normalized by library size, and all other default parameters were used. Analysis was performed for each cell type and stimulation condition. Evidence of association was quantified with the false discovery rate routine in DESeq2, and Q values of <0.05 were considered significant.

CD4+ T cell differentiation

CD4+ naïve cells were isolated from human PBMCs and murine splenocytes using species-specific Naïve CD4+ T Cell Isolation Kit (Miltenyi) and cultured for 5 days under polarizing conditions. Human cells were cultured with Dynabeads Human T-Activator CD3/CD28 (Invitrogen) with IL-2 (10 U/ml; PeproTech) supplemented with the following for TH1, TH2, and TH17 polarizations, respectively: IL-12 (1 ng/ml; PeproTech) and anti–IL-4 (10 mg/ml; clone MP4-25D2, BioLegend); IL-4 (4 μg/ml; PeproTech), anti–IFN-γ (1 μg/ml; clone B27, BioLegend), and anti–IL-12p70 (clone 20C2, eBioscience); and IL-1β (10 ng/ml; R&D Systems), IL-6 (10 ng/ml; PeproTech), IL-23 (10 ng/ml; R&D Systems), transforming growth factor–β1 (TGF-β1) (1 ng/ml; R&D Systems), anti–IL-4 (1 μg/ml), and anti–IFN-γ (1 μg/ml). Murine cells were cultured with anti-CD3 (5 μg/ml; clone 145-2C11) and soluble anti-CD28 (37.51) (1 μg/ml; eBioscience) supplemented with the following for TH1, TH2, and TH17 polarizations, respectively: IL-12 (1 ng/ml; PeproTech) and anti–IL-4 (10 mg/ml; clone 11B11, eBioscience); IL-4 (4 μg/ml; PeproTech) and anti–IFN-γ (10 μg/ml; clone XMG1.2, eBioscience); and IL-6 (20 ng/ml; PeproTech), IL-23 (10 ng/ml; R&D Systems), TGF-β1 (1 ng/ml; R&D Systems), and anti–IL-4 (10 mg/ml). For intracellular cytokine detection, cells were plated in fresh medium, stimulated with phorbol 12-myristate 13-acetate (PMA) (10 ng/ml) and ionomycin (1 ng/ml), and treated with brefeldin A (10 ng/ml; Sigma-Aldrich) and BD GolgiStop (BD Biosciences) for 5 hours before staining.

Enzyme-linked immunosorbent assays

Cell culture supernatants were assayed for human or murine IFN-γ, IL-4, IL-17A/F, and murine GM-CSF (granulocyte-macrophage colony-stimulating factor) using LEGEND MAX ELISA kits (BioLegend). Human plasma IgE was captured using anti-human IgE (clone G7-18, BD Biosciences). Human IgE full-length protein (Abcam) was the standard, and biotinylated anti-human IgE (clone G7-26, BD Biosciences) and Europium-conjugated streptavidin (PerkinElmer) were used for detection. Murine plasma IgE was captured using anti-mouse IgE (clone R35-72), with purified mouse IgE as the standard and biotinylated anti-mouse IgE (clone R35-118) for detection (BD Biosciences).

Animals

Tyk2 Ala1124–encoding targeting construct was generated by the Gene Recombineering Facility (Monash University). Gene targeting was performed in the embryonic stem (ES) cell line JM8F6. Successful homologous recombination between vector and endogenous Tyk2 locus was detected by long-range PCR using the 5HR-F2 (5′-GGGGCTCCTGTCTACAGCTC-3′) and 5HR-R2 (5′-TGTCGATCAGGATGATCTGG-3′), and the 3HR-F1 (5′-CGCGGGGATCTCATGCTGGAGTTC-3′) and 3HR-R1 (5′-CCCATGTGTGTCCCTTGATCCTGC-3′) primers and the LongAmp Taq PCR kit (New England Biolabs). Positive clones were further validated by Southern blotting using standard methods with probe amplification using the 5′-GGATCCGAACATGGTCCCTTGGATGT-3′ and 5′-GCTAGCGGCTTCACCTGTATCCCACT-3′ primers. Correctly targeted ES cells were injected into albino C57BL6/J blastocysts. Confirmation of the presence of transgenic Tyk2 at the transcript level was performed by PCR and Sanger sequencing using the 5′-CTGCTGGACAACGACAGG-3′ and 5′-ACACGCTGAACACGGAA-3′ primers.

Real-time quantitative PCR

Murine and human glyceraldehyde-3-phosphate dehydrogenase (GAPDH) housekeeping gene and murine Tyk2 and human IFNAR1 transcripts were detected using TaqMan Gene Expression Assays Mm99999915_g1, Hs02758991_g1, Mm00444469_m1, and Hs01066118_m1 (Applied Biosystems). Relative transcript levels are expressed as 2−∆Ct, where ∆Ct = (cycle threshold for the transcript of interest) – (cycle threshold for the housekeeping gene transcript).

The EAE mouse model

We immunized 8- to 12-week-old female mice subcutaneously in the flanks with 200 μg of MOG35–55 (Schafer) emulsified in incomplete Freund’s adjuvant (Sigma-Aldrich) supplemented with Mycobacterium tuberculosis H37Ra (4 mg/ml) (Difco; CFA) or with the pre-prepared Hooke MOG35–55/CFA immunisation kit for all cellular analyses (Hooke Laboratories). We administered 200 ng of Bordetella pertussis toxin (Sigma-Aldrich or Hooke Laboratories) intravenously on the day of immunisation and 48 hours later. The following EAE scoring system was used: score 0, no overt disease; score 1, weak tail; score 2, hindlimb weakness, abnormal gait, and/or impaired righting reflex; score 3, partial hindlimb paralysis; and score 4, complete hindlimb paralysis.

Murine cell isolation and stimulation

Murine peritoneal cells were plated in Dulbecco’s modified Eagle’s medium/F-12 with GlutaMAX (Gibco) and supplemented with 10% fetal calf serum (FCS); adherent cells were cultured for 24 hours before stimulation with IFN-β (10 ng/ml; Pestka Biomedical Laboratories), IFN-γ (10 ng/ml; PeproTech), or both for 18 hours. CNS tissue obtained from perfused mice was digested in 1× Hanks’ balanced salt solution (Sigma-Aldrich) supplemented with 0.05% collagenase D (Roche), tosyl-l-lysyl-chloromethane hydrochloride (TLCK) (0.1 μg/ml; Sigma-Aldrich), and deoxyribonuclease I (0.025 U/ml; Thermo Fisher Scientific). Mononuclear cells were then isolated by discontinuous Percoll (GE Healthcare Life Sciences) gradient centrifugation. Spleens and inguinal lymph nodes were obtained by standard methods. For culturing, cells were plated in RPMI supplemented with 10% FCS, 2 mM glutamine, penicillin/streptomycin, and with 1 mM sodium pyruvate, 0.1 mM nonessential amino acids, and 0.5 μM β-mercaptoethanol (Gibco). For phosflow, splenocytes were stimulated as for human PBMCs but with murine IFN-β (Pestka Biomedical Laboratories), IL-6, IL-10, IL-12 (PeproTech), or IL-23 (R&D Systems). For recall assays, cells were stimulated for 72 hours with MOG35–55 (10 μg/ml) and restimulated with PMA and ionomycin as for CD4+ T cell differentiation.

Statistical analysis of HEK cells, human immune cells, and murine data

Statistical tests were performed using GraphPad Prism and the R statistical software package. For analyses comparing two unpaired experimental groups, Mann-Whitney U test was used, whereas for paired analyses, Wilcoxon matched-pairs rank test was used, with a 5% significance level, and estimating one- or two-tailed P values depending on previous assumptions regarding the directionality of the effect tested. For unpaired primary human pSTAT analyses, Bonferroni correction for multiple testing was used (taking P = 0.05 and considering eight independent hypotheses for four different cytokine types and two SNPs), and the significance threshold estimated was P = 0.00625, noting that different cell subsets were not considered as independent tests because all are TYK2-positive and the SNPs are missense such that different subsets were analyzed to account for differences in cytokine receptor levels that might confound analyses of bulk populations. At this significance threshold, we obtained >80% power to detect the differences observed. For these analyses by rs34536443 and rs12720356 genotype, P values provided are for nonadditive and additive genetic models, respectively, noting that age and sex were not included as covariates because no age or sex effects were observed (all P > 0.05).

Structural analysis

Site-directed mutagenesis was applied to a previously reported JAK2 kinase domain construct (50) to produce a version bearing P1057A. DNA encoding the Ala1057 JAK2 kinase domain was manufactured by GenScript. Protein was expressed and purified as previously described (50). Purified Ala1057 JAK2/TG101209 complex [6 mg/ml in 250 mM NaCl, 25 mM bicine (pH 8.5), 0.5 mM tris(2-carboxyethyl)phosphine (TCEP), and 1 mM inhibitor (from 50 mM dimethyl sulfoxide)] was subjected to crystallization trials using screening solutions from Hampton Research and Qiagen. Optimized conditions used a reservoir with 1.9 to 2.0 M sodium malonate (pH 7.0), 0.1 M bis-tris propane (pH 7.0), and 10 mM ZnCl2 in hanging drops made from 2 μl of protein solution plus 2 μl of reservoir at 18°C. Crystals were harvested and placed briefly in 100% reservoir as cryoprotectant and then plunged into liquid nitrogen. Diffraction data were collected at 110 K using APS beamline 21-ID-G (wavelength, 0.97856 Å). Data were integrated and scaled using HKL2000 and elements of the CCP4 suite. Initial phases were obtained by molecular replacement with Phaser using PDB 4JI9 (JAK2) as search probe. Refinement was performed using Refmac5, phenix.refine, and BUSTER-TNT (final model, 94% favored Ramachandran angles), and model manipulation and electron density inspection used Coot.

U.K. Biobank data analysis

U.K. Biobank genotype data included 152,729 subjects. We excluded individuals of non-European ancestry (based on principal components analysis) and samples with abnormal heterozygosity rates and high missing genotype rates, leaving 116,732 individuals [male/female ratio of 61,339:55,214; mean age (±SEM) = 63.82 (±7.94) years]. Data on hospital episode statistics contained in-patient data from 01 April 1997 onward and included ICD-10–coded data from admissions, operations, and procedures. We found no association between genotype and number or nature of hospitalization events, aside from associations with autoimmune disease. In TYK2-deficient individuals, the effect on immunodeficiency is 100% penetrant (29). In the U.K. Biobank, 1.2% of the cohort had at least one hospital episode recorded with relevant infections. With this prevalence and with 249 rs34536443 C/C individuals, we had >80% power to detect an association with a relative risk of >6.0 (using a recessive genetic model with a minor allele frequency of ~4%) and a case-control ration of 10 (at α = 0.05). From the cancer registry, genetic data were available for 13,315 subjects, and the most common cancer type (melanoma and other malignant neoplasms of the skin) has a prevalence of 3.4%. With this prevalence, we had >80% power to detect an association (using a recessive genetic model as above) with a risk ratio of 2.621. Survival analysis (using the Cox proportional hazard regression model) with rs34536443 genotypes, and correcting for the effect of gender in life expectancy, resulted in no significant differences. In the U.K. Biobank, there were 2146 death events registered, and we estimated power to detect a genetic association (recessive model) with the survSNP R package, assuming parameters obtained from the U.K. Biobank data. The study was powered to detect an association with reduced survival, with a genotype hazard ratio of 2.6 (considerably above the effect estimated for gender, 1.5). All data used were last updated June 2015.

SUPPLEMENTARY MATERIALS

www.sciencetranslationalmedicine.org/cgi/content/full/8/363/363ra149/DC1

Fig. S1. CRISPR-Cas9–modified HEK293T cell lines and their respective TYK2 expression and phosphorylation.

Fig. S2. CRISPR-Cas9–modified TYK2 knockout cell lines, their respective TYK2 expression, and impact of genetic modification on cytokine signaling.

Fig. S3. Summary data for Oxford BioBank blood donor cohort and representative optimization data for cytokine stimulation.

Fig. S4. Representative cytokine-induced pSTAT and IFN-β–induced pSTAT1 by the rs34536443 genotype in human primary immune cells.

Fig. S5. Downstream consequences of the rs34536443 genotype–dependent differential cytokine signaling.

Fig. S6. Signal 3 haplotype colocalization with epigenetic marks.

Fig. S7. Characterization of Tyk2 Ala1124 mice.

Fig. S8. Tyk2 P1124A mouse weight change and incidence after immunization, and draining lymph node cell cytokine production after restimulation with MOG peptide.

Fig. S9. Cytokine production after MOG peptide restimulation of CNS CD4+ T cells and clinical profile of Tyk2 Pro/Ala1124 mice after immunization.

Fig. S10. Immune profile of Oxford BioBank donor pairs by the rs34536443 genotype.

Fig. S11. Cytokine receptor expression by the rs34536443 genotype.

Fig. S12. Alignment of human TYK2 and JAK2 kinase domains and statistics for JAK2 (P1057A)/TG101209 structure.

Table S1. Summary of previously reported genetic associations in the TYK2 region.

Table S2. TYK2 region 90% credible sets for the three independent association signals identified by multinomial association analysis.

Table S3. TYK2 region associations for each of the autoimmune diseases included in the multinomial association analysis.

Table S4. Summary of TYK2 transcript expression by rs34536443 and rs9797854 genotypes across immune cell types.

Table S5. Evidence for nonadditivity of rs34536443 effect on pSTAT.

Table S6. Analysis of normalized pSTAT mean fluorescence intensity by the rs34536443 genotype.

Table S7. Analysis of frequency of pSTAT-positive primary human immune cells by the rs12720356 genotype.

Table S8. Analysis of normalized pSTAT mean fluorescence intensity in primary human immune cells by the rs12720356 genotype.

Table S9. Summary of eQTL colocalizations with TYK2 region association signal 3.

Table S10. Association of the rs34536443 genotype with U.K. Biobank health record phenotypes.

Table S11. Relationship between genotype-dependent cytokine signaling and risk for autoimmune disease and certain recurrent or severe infections.

REFERENCES AND NOTES

  1. Acknowledgments: We thank the Genetic Analysis of Psoriasis Consortium, the International Genetics of Ankylosing Spondylitis Consortium, the International Multiple Sclerosis Genetics Consortium, the U.K. IBD Genetics Consortium, and the Wellcome Trust Case Control Consortium 2 for data access. We thank the volunteers from the Oxford BioBank, National Institute for Health Research (NIHR) Oxford Biomedical Research Centre (BRC), for their participation. The Oxford BioBank (www.oxfordbiobank.org.uk) is also part of the NIHR National Bioresource, which supported the recalling process of the volunteers. We also thank A. Smith and T. Sauka-Spengler for help in establishing the Southern blotting and ATAC-Seq methods, respectively; C. Monaco, P. Amjadi, and D. Ahern (Kennedy Institute of Rheumatology, University of Oxford) for mass cytometer access and technical support; and B. Davies and D. Biggs (Transgenic Core, Wellcome Trust Centre for Human Genetics, University of Oxford) for transgenic mouse generation services. This study makes use of data generated by the Blueprint Consortium; a full list of the investigators who contributed to the generation of the data is available at www.blueprint-epigenome.eu, and funding for the project was provided by the European Union's Seventh Framework Programme (FP7/2007-2013) under grant agreement number 282510 BLUEPRINT. This research has also been conducted using the U.K. Biobank Resource. We would also like to acknowledge the flow cytometry facility at the Weatherall Institute of Molecular Medicine (WIMM), which is supported by the Medical Research Council (MRC) Human Immunology Unit, the MRC Molecular Haematology Unit (MC_UU_12009), the NIHR Oxford BRC and John Fell Fund (131/030 and 101/517), the EPA fund (CF182 and CF170), and the WIMM Strategic Alliance awards (G0902418 and MC_UU_12025). We thank the High-Throughput Genomics Group at the Wellcome Trust Centre for Human Genetics (funded by Wellcome Trust grant reference 090532/Z/09/Z) for the generation of sequencing data. We thank Shamrock Structures, LLC for diffraction data collection. Use of the APS, an Office of Science User Facility operated for the U.S. Department of Energy (DOE), Office of Science by Argonne National Laboratory, was supported by the DOE under contract no. DE-AC02-06CH11357. Use of the LS-CAT Sector 21 was supported by the Michigan Economic Development Corporation and the Michigan Technology Tri-Corridor (grant 085P1000817). Funding: Work in the authors’ laboratories is supported by the U.K. and Danish MRCs, the Alan and Babette Sainsbury Charitable Fund, the Naomi Bramson Trust, the Clinical Neuroimmunology Fund, the Oxford BRC, the Oak Foundation (L.F.), the Rosetrees Trust (L.F. and C.A.D.), and the Wellcome Trust (100308/Z/12/Z and 106130/Z/14/Z to L.F. and 100956/Z/13/Z to G.M.). L.S. was supported by a Christopher Welch Scholarship. Author contributions: C.A.D. contributed to the conception, coordination, and design of the study, all experiments apart from the structural biology, and drafting and writing of the manuscript. A. Cortes performed all statistical genetics, epigenetic and next-generation sequencing, and U.K. Biobank data analyses and contributed to the writing of the manuscript. L.S. contributed to the blood sample processing and analyses and the study conception and design, and L.S. and S.L.F. contributed to preliminary immunoblot analyses. H.G.E. generated the CRISPR-Cas9–modified cell lines, contributed to blood sample processing, and performed the FACS sorting. K.E.A. and G.K. performed the animal immunizations, K.E.A. contributed to murine cell activation experiments and performed all murine perfusions, and C.D. contributed to flow cytometry of murine samples. L.J. contributed to the design and execution of genetic and statistical analyses. T.B. and O.A.L. contributed to blood sample processing. S.B.K. contributed to TaqMan SNP genotyping. J.C., M.J.N., and F.K. contributed to the coordination and sample obtainment from the Oxford BioBank, and M.J.N. performed genotyping for volunteer preselection. S.S. and A. Compston contributed multiple sclerosis genetic and clinical data. M.U., A.R.J., C. Everett, and C. Eigenbrot designed and performed the crystallography and relevant biophysical analysis and contributed to the writing of the final manuscript. J.I.B. contributed to the writing of the final manuscript. G.M. contributed to the conception and design of genetic and statistical analyses and to the writing of the manuscript. L.F. contributed to conception, design, and coordination of the study and drafting and writing of the manuscript. Competing interests: L.J. is a consultant to Genomics plc. M.U., A.R.J., C. Everett, and C. Eigenbrot hold equities in Roche Holdings. J.I.B. is on the board of Roche and holds equities in Roche Holdings. G.M. is a cofounder of, holder of shares in, and consultant to Genomics PLC. The other authors declare no competing financial interests. Data and materials availability: The x-ray crystal structure of the Ala1057 JAK2 kinase domain at 2.66 Å resolution has been deposited in the PDB (PDB accession 5HEZ).
View Abstract

Navigate This Article