Natural Killer cells demonstrate distinct eQTL and transcriptome-wide disease associations, highlighting their role in autoimmunity

Abstract

Natural Killer cells are innate lymphocytes with central roles in immunosurveillance and are implicated in autoimmune pathogenesis. The degree to which regulatory variants affect Natural Killer cell gene expression is poorly understood. Here we perform expression quantitative trait locus mapping of negatively selected Natural Killer cells from a population of healthy Europeans (n = 245). We find a significant subset of genes demonstrate expression quantitative trait loci specific to Natural Killer cells and these are highly informative of human disease, in particular autoimmunity. A Natural Killer cell transcriptome-wide association study across five common autoimmune diseases identifies further novel associations at 27 genes. In addition to these cis observations, we find novel master-regulatory regions impacting expression of trans gene networks at regions including 19q13.4, the Killer cell Immunoglobulin-like Receptor region, GNLY, MC1R and UVSSA. Our findings provide new insights into the unique biology of Natural Killer cells, demonstrating markedly different expression quantitative trait loci from other immune cells, with implications for disease mechanisms.

Fig. 1. cis eQTL in NK cells.

a Total significant (FDR < 0.05) primary and conditional cis eQTL in NK cells from 245 individuals. Numbers of eQTL with evidence of colocalisation with at least one GWAS trait (RTC > 0.9) are highlighted (pink). b Frequencies of cis eQTL (FDR < 0.05) specific to NK cells or shared with other immune cells; monocytes, neutrophils, CD4⁺ and CD8⁺ T cells. c Enrichment of shared causal loci between NK cell eQTL (total, left panel; NK cell-specific, right panel) and GWAS traits (n = 100). Traits are compared with the enrichment observed for height as background. Significant traits (FDR < 0.05) are highlighted (pink). Point size is proportional to a trait's number of GWAS-significant loci. d Effect of rs1788098 genotype on CD226 expression in NK cells (n = 245 individuals). Box and whisker plot; boxes depict the upper and lower quartiles of the data, and whiskers depict the range of the data excluding outliers (outliers are defined as data-points > 1.5 × the inter-quartile range from the upper or lower quartiles). e The effect of rs1788098 genotype on CD226 expression is specific to NK cells. Significant eQTL effects are highlighted (green). f The CD226 eQTL in NK cells colocalises with a risk locus for type-1 diabetes. SNPs are coloured according to strength of LD (CEU population) to the peak eSNP (rs1788098); brown r² > 0.8, orange 0.5 < r² ≤ 0.8, yellow 0.2 < r² ≤ 0.5, grey r² ≤ 0.2. g Association of rs1788097 (exact proxy for rs1788098 in European populations, r² = 1) with autoimmune diseases and haematological indices. GWAS-significant (p < 5 × 10⁻⁸) associations are highlighted (pink). h Regional association plot of the primary IRF5 eQTL in NK cells. i The primary IRF5 eQTL does not colocalise with a GWAS risk locus in the IRF5 region for systemic lupus erythematosus. j Conditioning on the peak IRF5 eSNP reveals an independent, secondary eQTL for IRF5 in NK cells. k The secondary IRF5 eQTL colocalises with the IRF5 region for systemic lupus erythematosus. Genotype to phenotype correlations were calculated with linear regression. P values are two-sided.

Fig. 2. trans eQTL in NK cells.

a Manhattan plot depicting significant trans eQTL in NK cells. Physical position (x axis) represents location of the target gene. Coloured points represent significant (FDR < 0.05) trans eQTL; 2266 trans SNP-gene associations, affecting expression of 64 genes. Regulatory networks are highlighted as follows; GNLY (blue), MC1R (orange), UVSSA (yellow), KIR (pink), all others (green). b Effect of rs1866140 genotype on GNLY expression in NK cells (n = 245 individuals). c A trans-regulatory network of 5 genes, mediated in cis by GNLY expression. d Effect of rs1866140 genotype on GNLY regulatory network genes in trans (n = 245 individuals). PP4, posterior probability of a shared causal variant with GNLY cis eQTL. e Effect of rs117406136 genotype on MC1R expression in NK cells (n = 245 individuals). f A trans-regulatory network of 4 genes, mediated in cis by MC1R expression. g Effect of rs117406136 genotype on MC1R regulatory network genes in trans (n = 245 individuals). PP4, posterior probability of a shared causal variant with MC1Rcis eQTL. h Effect of rs111632154 genotype on UVSSA expression in NK cells (n = 245 individuals). i A trans-regulatory network of 4 genes, mediated in cis by UVSSA expression. j Effect of rs111632154 genotype on UVSSA regulatory network genes in trans (n = 245 individuals). PP4, posterior probability of a shared causal variant with UVSSA cis eQTL. k A trans-regulatory network of 10 genes, mediated mediated by KIR2DS4del. l Effect of KIR2DS4del genotype on KIR regulatory network genes in trans (n = 245 individuals). Box and whisker plots; boxes depict the upper and lower quartiles of the data, and whiskers depict the range of the data excluding outliers (outliers are defined as data-points > 1.5 × the inter-quartile range from the upper or lower quartiles). Genotype to phenotype correlations were calculated with linear regression. P values are two-sided.

Fig. 3. NK cell TWAS in autoimmune diseases.

Manhattan plot of autoimmune disease TWAS. All significant trait-associated genes (n = 98) are labelled. Novel gene-trait associations (no GWAS-significant locus within 1Mb and no evidence supporting gene-trait association in Open Targets Genetics) are highlighted (yellow, n = 27). Novel gene-trait associations which colocalise (posterior probability colocalisation > 0.8) with an NK cell-specific eQTL are circled (pink, n = 3).

Fig. 4. NK cell eQTL determine protein expression on NK cells in patients with metastatic melanoma.

Genotype-phenotype correlations are calculated by linear regression. P values are calculated with likelihood ratio tests. Expression of CD57 (orange, n = 78) and CD226/DNAM1 (blue, n = 24) protein on NK cells is determined by the respective RNA eSNPs in patients with metastatic melanoma. KIR2D antigen expression on NK cells in patients with metastatic melanoma (n = 78) is modified by KIR2DS4del copy number (pink). Box and whisker plots; boxes depict the upper and lower quartiles of the data, and whiskers depict the range of the data excluding outliers (outliers are defined as data-points > 1.5 × the inter-quartile range from the upper or lower quartiles). Genotype to phenotype correlations were calculated with linear regression. P values are two-sided.