|Home | About | Journals | Submit | Contact Us | Français|
Natural killer (NK) cells are now recognized to exhibit characteristics akin to cells of the adaptive immune system. The generation of adaptive memory is linked to epigenetic reprogramming including alterations in DNA methylation. The study herein found reproducible genome wide DNA methylation changes associated with human NK cell activation. Activation led predominately to CpG hypomethylation (81% of significant loci). Bioinformatics analysis confirmed that non-coding and gene-associated differentially methylated sites (DMS) are enriched for immune related functions (i.e., immune cell activation). Known DNA methylation-regulated immune loci were also identified in activated NK cells (e.g., TNFA, LTA, IL13, CSF2). Twenty-one loci were designated high priority and further investigated as potential markers of NK activation. BHLHE40 was identified as a viable candidate for which a droplet digital PCR assay for demethylation was developed. The assay revealed high demethylation in activated NK cells and low demethylation in naïve NK, T- and B-cells. We conclude the NK cell methylome is plastic with potential for remodeling. The differentially methylated region signature of activated NKs revealed similarities with T cell activation, but also provided unique biomarker candidates of NK activation, which could be useful in epigenome-wide association studies to interrogate the role of NK subtypes in global methylation changes associated with exposures and/or disease states.
Epigenetic modifications, including DNA methylation, dictate programmed lineage differentiation within the immune system.1 Remodeling the epigenome during development leads to progressively restricted immune subtypes. DNA methylation is thought to provide a chemically stable mark that serves to initiate and reinforce these cell fate decisions.2 Large-scale methylation studies of specific lineage pathways have been performed on differentiating erythroid,3,4 T lymphocyte,5,6 B lymphocyte,7 and myeloid4 cells. These important developmental studies reveal insights into the regulation of cell fate decisions but also serve as a potential source of biomarkers that help characterize different immune compartments during the immune response. The best known example of a single gene epigenetic marker of immune cell lineage is demethylation within the FOXP3 gene that provides a marker of CD4 T regulatory cells (Treg).8 FOXP3 expression drives the transcriptional program of Tregs.5 Another example of lineage specific methylation, which marks human natural killer cells (NK), is represented by demethylated CpG sites within the NK cell activation receptor NCR1/NKp46 gene locus that predicts the numbers and proportions of NK cells in blood and tissue.9
Less understood is the extent to which stable and heritable changes in DNA methylation might be produced during immune activation within cells of a common lineage identity. Although transcriptomic studies implicate genome wide alterations in gene expression during pathogen- and cytokine-induced immune activation in different cell types,10-14 the exact role of DNA methylation in the genesis of these effects is not known. However, a few individual loci have been studied intensively. For example, the production of interferon γ (IFNG) in response to IL2 treatments or T cell receptor (TCR) activation leads to the demethylation of the IFNG locus in CD4 Th1 cells.15,16 Studies of CD8 T cell activation and memory generation point to a profound plasticity in DNA methylation within the CD8A TCR co-receptor locus,17 which correlates with different cytokine exposures and conditions of Th1/2 polarization. Studies using broader CpG scans showed that CD4 helper T cells undergo dynamic and flexible methylation changes both in a genome-wide and a targeted manner during T cell activation and differentiation.18 In addition, ex vivo differentiation of human CD14 blood monocytes to immature and mature dendritic cells led to significant remodeling of the monocyte/macrophage methylome.19
Although the NK cell methylome has not been extensively studied during cytokine or NCR receptor activation, some work has been done with cytomegalovirus (CMV) infected mice and humans. A clonal expansion of NKG2C receptor positive NK cells has been shown to occur, and in some HCMV infected individuals, subsets of memory-like NK cells deficient in B-cell and myeloid signaling proteins, including the tyrosine kinase SYK, have been identified. Furthermore, the downregulation of SYK occurred in the context of DNA hypermethylation within a restricted segment of the SYK gene promoter.20 HCMV-associated NK cells also displayed deficiencies and hypermethylation in signaling adaptors such as EAT-2 and FCER1G as well as the ZBTB16/PLZF transcription factor.21
We predicted that remodeling of the NK methylome following NKp46 engagement would be a feature of NK activation based on the prior demonstration that NK cells display a type of immunologic memory that shares many features with CD8 T effector and memory cells following antigen activation.22 The work herein studied the DNA methylation profiles of resting and activated NK cells; NK activation was achieved by ex vivo antibody-ligation of the CD2 and NKp46 receptors and cell culture in the presence of rIL2. We note similarities in the induced methylation alterations with those occurring in other antigen and cytokine stimulated cell lineages.
A specific aim of the study was to identify promising single gene candidates that might provide biomarkers of NK cell activation. We developed a single gene qPCR assay of demethylation within the BHLHE40 transcription factor using the highly sensitive droplet digital PCR (ddPCR) technique23 and consider its potential as a biomarker of NK activation within peripheral whole blood.
Untouched human NK cells were isolated from a total of 12 individuals over 2 phases of this study (Suppl Fig. 1). Isolated naïve and activated NK cell populations were confirmed by flow cytometric analyses (Suppl Fig. 2). On average, 96.1% (±2.7%) of NK cells were activated at day 9 (n = 12; data not shown) using CD54 as a measure of activation, compared with no activation at T = 0. The killing phenotype of activated cells was confirmed using a target cell lysis assay and the ability of activated cells to demonstrate antibody-dependent cell cytotoxicity measured in the presence and absence of the monoclonal antibody therapeutic agent Rituximab (Fig. 1).
The experimental design of the study was based on a discovery phase followed by a replication phase (Suppl Fig. 1). During the discovery phase, isolated human NK cells were harvested at 3 time points: just after MACS isolation (T = 0), 16 h after initiation of activation (T = 16 h), and 9 d post-activation (T = 9 d). Cell viability, purity, and activation status were verified using flow cytometric techniques (Suppl Fig. 2 and Suppl Table 1). Hierarchical clustering analysis of the discovery phase array data indicated there was no significant difference between the methylation patterns at T = 0 and 16 h (data not shown). Somewhat surprisingly, this included no indication of any altered DNA methylation status for genes known to be involved in cell cycle regulation. Hence, the replication studies focused on T=0 and 9 d post-activation.
NK cells were isolated from 6 additional donors for the replication phase of the study. As with the discovery set, these were analyzed by hierarchical clustering, which was able to determine differential methylation patterns between naïve and activated NK cells (Fig. 2). Additionally, the genomic location of the top 1000 loci was conserved across the discovery and replication data sets (data not shown). The discovery and replication datasets were compared in order to generate a list of the top significant loci in common between these 2 sets (Suppl Table 3). This analysis returned 385 loci containing 259 identified genes, microRNA, and coding regions. Fig. 3 shows the genomic distribution characteristics and locations of the significant probes within gene regions.
In an effort to more completely investigate the nature of the altered loci observed in NK activation, we interrogated the lineage of origin of the top 500 replicated loci by applying the algorithm of Houseman24 to a well-known library of the DNA methylation profiles of the common cell types in peripheral blood.25 This algorithm has been validated as a tool to discern lineage of unknown peripheral blood mixtures, using known profiles from isolated cells of a single lineage.26 This analysis revealed that there were 7 DMSs that were found in activated NK cells that were not differentially methylated in other cell types, including cg11783497 (IL1RN), cg07996532 (ECE1), cg17566874 (CSF2), cg05550420 (DCHS1), cg06169746, cg26174398, cg13157980 (EIF2C2). This implies some level of specificity of DNA methylation in activated NK cells (although these genes may be expressed in other immune cell types).
To identify the most phenotypically important loci, a series of prioritizing criteria were developed (Suppl Fig. 1B) that weighted effect size and regional coherence (multiple significant probes). This identified several genes with considerable differential methylation between naïve and activated NKs; these genes were further analyzed for altered gene expression. These included genes that: i) ranked highest on the priority list (including BHLHE40, MYO1G, NCF4); ii) were reported by others to have altered expression in NKs treated with CMV (such as BHLHE40, KSR1, MAML2, TNFAIP3, ZBTB32);11 iii) appear on our list of top 27 genes (e.g. BHLHE40, IL13, MYO1G and NCF4); and/or iv) genes known to play a role in activation of other immune cells (BHLHE40, KSR1, NCF4). All of the candidates tested displayed significant hypomethylation and significantly increased mRNA expression was observed for 6 of the tested genes while a small but significant decrease in expression was found for one. The remaining 3 genes exhibited no change in expression (Suppl Fig. 4).
Bioinformatics analyses of the 285 most significant coding sequence-associated loci in common between the discovery and replication data sets were assessed for functional annotation enrichment and for pathway associations. Of these 285 coding sequence-associated loci, 278 were successfully mapped by DAVID and generated 15 significant, functionally similar clusters (enrichment score > 1), with the highest enrichment for regulation of leukocyte/lymphocyte/cell activation (enrichment score of 4.03; Table 1A). The functional annotation chart module returned 164 significant (P < 0.05) overrepresented terms (Table 1B). Approximately 50% of these are directly immune function related (e.g., cytokine/cytokine receptor interaction, regulation of leukocyte activation, immune response) while many of the remaining are clearly basic cellular functions, but likely require enhanced activity during cell activation (e.g., regulation of phosphorylation, protein transport and protein secretion, etc.) Additionally, Ingenuity Pathway Analysis using the “canonical pathway” module revealed significant enrichment of immune-related signaling pathways (e.g., NF-kappa β, TNFR2 and IL-10 signaling pathways and cytokine production pathways), differentiation pathways (e.g., T-helper and dendritic cell maturation), as well as disease-related pathways with key immune cell components (e.g., asthma and multiple sclerosis). Significant (Benjamini-Hochberg P-value < 0.05) immune related pathways were observed as the top hits regardless of whether the analysis was constrained by immune cell filters or not (data not shown).
The frequency of differential methylation within the transcription factor binding sites for 161 transcription factors was assessed using bioinformatics techniques to mine the ENCODE database (Suppl Table 6). We hypothesized there would be enrichment for a limited number of transcription factors that might suggest their importance in mediating the activation phenotype. Separate exploration of significant hypo- and hypermethylated DMS revealed a surprisingly similar list of transcription factor binding sites perturbed by changes in methylation. Significant associations involved EP300, CEBPB, FOS, MYC, GATA2, STAT3, RUNX3 and others. EP300, STAT3, and RUNX3 play well-known roles in NK differentiation; however, they are also expressed in other immune cell lineages. A novel observation in our data was the common occurrence of GATA2 among the TFBSs; GATA2 has recently been shown essential in the terminal maturation of CD56 dim NK cells that are most enabled for cytotoxicity (Suppl Table 6). Two TFBS perturbations were unique to the hypermethylated CpG: GRp20 (glucocorticoid receptor) and ZKSCAN1. There were 16 TFBS that had hypomethylated DMS that were not altered by hypermethylation.
In other functional analyses, the 100 non-coding loci were evaluated using the GREAT software application (Suppl Table 7). Results of these analyses align with alterations in immune activity, including iNKT development, T cell effector and memory function and STAT3 signaling. This analysis also revealed enrichment of targets of miR-155, including SOCS1. Recent studies show that miR-155 is upregulated in NK cells by activating receptor ligation and cytokines (IL-12, IL-15, IL-18) and regulates INFG production. Interestingly miR-155 is involved in the expansion of virus-specific NK cells, the latter being mediated through specific inhibition of SOCS1.27 Taken together, these annotation analyses confirm the importance of the differentially methylated loci identified in our study as playing an important role in NK cell activation.
Further analysis was made by comparison of the most significant 385 differentially methylated loci associated from this study to a known activated NK model cell line, the genetically altered NK-92 cells (Conkwest). Unlike the parental NK-92 strain, the CD16.176V.NK-92 line has been transfected to express the high affinity variant of the CD16 Fc receptor, making them optimized for ADCC. A hierarchical clustering analysis of freshly isolated naïve NK, experimentally activated NK, and genetically modified NK-92 strain showed, as expected, the NK-92 to cluster more closely with the chemically activated NKs than naïve cells (Fig. 4). To better understand the similarity between these cells optimized for ADCC and the in vitro activated NK cells, a comparison of significant differentially methylated sites from both groups was performed. The analysis revealed an overlap of 169 loci, 111 of which were associated with gene/coding sequence and 78 loci (representing 58 coding sequences) unique to NK-92 (Suppl Table 4). Based on the DAVID functional enrichment analysis of this overlapping dataset, the most significant enriched functional annotations were biological processes of immune response, inflammatory response, positive regulation of lymphocyte activation and regulation of lymphocyte activation (P-value range: 9.07 × 10−7 to 3.65 × 10−5; data not shown). The functional annotation clustering identified the most important members were associated with positive regulation of lymphocyte, leukocyte and cell activation (3 separate “GOTERM_BP_FAT” terms), with an enrichment score of 4.34 (data not shown). Enrichment analysis of the 58 sequence-associated DMSs for NK-92 revealed a significant enrichment of genes for cell death/apoptosis (enrichment score 1.75) and transcription regulation (enrichment score 1.22; data not shown).
A number of potential markers of activated NK cells were tested for suitability for a digital PCR assay. In an effort to rank order the changes in DNA methylation, the Illumina HumanMethylation450 BeadChip DMS candidates were prioritized by statistical significance (P- and q-values), extent of coverage (multiple positive probes), effect size, known biological role, and novelty (data not shown). Of the trial candidates, BHLHE40 showed the greatest potential for assay design. A region of exon 5 that included 3 CpGs interrogated by the array and 25 other CpGs, was bisulfite sequenced (Fig. 5). A combination of differential methylation and sequence amenable to assay design resulted in the use of 2 particular CpGs as the target region to assess the activation phenotype.
The newly developed BHLHE40 demethylation ddPCR assay was applied to the naïve and activated NK cells from 3 donors (of the 12 used in the NK study; Fig. 6). Additionally, 8 different isolated cell types were assessed to determine the specificity of BHLHE40 demethylation as a marker for activated NK cells (not subjects from this study; Fig. 6). The in vitro activated NK cells consistently measured ~60% demethylation as compared with ~15–35% demethylation from naïve NKs. Some demethylation was also detectable in other isolated cell types, primarily from T- and B-cells. The robustness of the assay was also assessed by measuring BHLHE40 demethylation in whole blood samples from 6 individual donors (not subjects from this study; Fig. 6 inset). The assay demonstrated linear kinetics in a 3 log dilution series (Rsquare > 0.98; data not shown).
The Illumina HumanMethylation450 BeadChip arrays were applied to paired donor primary naïve and in vitro activated NK cells in an effort to characterize the profile of DNA alterations associated with NK activation. This same methodology can be utilized to investigate and develop additional biomarkers of NK activation.
A significant observation from these experiments is that NK cells display no detectable changes in DNA methylation 16 h following activation with antibody mediated receptor ligation and IL2 stimulation. In addition, no DNA methylation changes in genes associated with cell cycle regulation were observed. The short-term kinetics of DNA methylation changes in stimulated immune cells has been controversial. While rapid changes in gene expression have been described in naïve T cells differentiating to effector and memory cells, little consensus exists on the methylation alterations that accompany these early events.10,13 Our observations at 16 h are consistent with previously published data studying CD4 T cells undergoing in vitro activation with CD3 and CD28 antibodies.6 Further studies are necessary to determine if the lack of early kinetic change in methylation is due to steps involving active TET-mediated hydroxymethylation or a passive dilution and lack of methylation maintenance. In contrast to the lack of short-term changes, NK cells at day 9 post-activation (this study) showed great numbers of highly significant changes in DNA methylation, reflecting a broad remodeling of the methylome.
A discovery and replication experimental design was employed, applying strict false discovery rate thresholds, that identified significant differential methylation (predominately hypomethylation) localizing to intragenic regions. We identified several hypomethylated loci that are among those previously shown to be differentially methylated in the immune system, whereas, many other loci discovered have never been associated with differential methylation in activated immune cells, including NKs. The functional annotation analyses confirmed the biological relevance of the differentially methylated loci to NK cell activation. The 111 gene/coding sequence-associated loci found in common between the in vitro activated human NK cells and the genetically modified NK-92 cell line enhanced for ADCC may represent candidate genes (and critical loci) that will be useful in pharmacologic manipulation to enhance NK activation.
A prominent focus of demethylation in NK cell activation involved the chromosome 6p21.3 region that contains the lymphotoxin α (LTA), tumor necrosis factor α (TNFA) and lymphotoxin β (LTB) genes. This 10 kb region has been intensively studied for epigenetic regulation of pivotal inflammatory mediators.28 Both LTA and TNFA were demethylated in our activated NK cells. Previous work showed that the DNA methylation status of the proximal promoter of the TNFA gene is highly correlated with gene transcription.28 Three sites within exon 1 of LTA in NK cells overlapped with loci demethylated within naïve CD4 T cells in Sjogren's syndrome.29
This study found additional members of the TNF family outside the 6p21 locus were also differentially methylated in activated NK cells. These include members of the TNF (ligand) superfamily, members 4, 8, and 10, and TNF receptor superfamily members 4 and 9. All of these exhibited significantly hypomethylated loci in activated NK cells except TNFSF8, which showed a significant hypermethylated DMS upon activation. Changes in methylation status for TNFRSF4, 9 and TNFA were observed during Treg differentiation, also exhibiting demethylation.30 Table 2 illustrates more examples of differentially methylated genes from other activated immune cells that appear to be common with NK cell activation.30-39
The Th2 cytokine locus has been intensively studied as a model of chromatin dynamics directing gene expression in the immune system. The genes encoding the Th2 cytokines IL4, IL5 and IL13 are clustered in a 160 kb region of chromosome 5q31 in humans. Significant methylation changes were observed in loci associated with IL5 and IL13 (Suppl Fig. 3). One of the strongest and most consistent differentially methylated gene loci in our experiments was the demethylation of CpG sites within 1 kb of the proximal promoter of the IL13 gene in activated NK cells. Eight consecutive probes displayed highly consistent demethylation in activated compared to naïve NK cells, which was associated with marked upregulation of mRNA synthesis (Suppl Fig. 4). Studies by others showed demethylation in this region was characteristic of human Th2 polarized CD4/CD28 activated T cells leading to increased IL13 synthesis.40
Our studies extend the association of immune activation with demethylation of the proximal IL13 promoter to NK cells. This is of interest as NK cells are typically thought of in relationship to Th1 inflammatory cytokines, such as IFNG. Upregulation of IFNG was confirmed in our NK cells by flow cytometry, however, NK cells also produce IL13 and other Th2 cytokines,41 a phenotype implicated in allergic disease.42 Additionally, there is evidence of the existence and regulation of 2 distinct resting peripheral NK cell subsets producing either Th1 or Th2 cytokines exclusively.43 Our studies are unable to discern whether DNA methylation states of the IFNG and IL13 genes are coordinated within individual NK cells.
Also located at 5q31 is CSF2, the gene coding for the cytokine that controls production, differentiation, and function of granulocytes and macrophages. We found 3 CpGs sites reproducibly demethylated at day 9 including a site (cg02325250) previously identified as being rapidly demethylated in PMA/I stimulated mouse CD4 cells and associated with enhanced CSF2 expression.18 However, in our experiments at 16 h post-stimulation, none of the probes in CSF2 were significantly different from those of naïve NK cells, suggesting different kinetics in NK and T cells for this gene locus. By 9 d, NK cells displayed consistently lower methylation. We note that BHLHE40 activity is required for CSF2 production, and both genes are targets of demethylation in activated NK cells.44
Demethylation of numerous genes encoding established (for other immune cell types) immune stimulatory and activation loci in NK cells underscores the common epigenetic programs that are shared between the innate and adaptive immune systems. This is particularly true with respect to NK and T cells,11 and is further supported by findings in this study. Table 2 provides several examples of stimulatory and activation-related genes in T cells that also exhibit significant differential methylation during NK cell activation.
The IFNG gene was broadly demethylated in resting NK cells with no statistically significant genome-wide differential methylation following activation, although a small increased IFNG expression was observed by FACs analysis. Our DNA methylation result is consistent with previous studies showing that IFNG production in peripheral NK cells does not require DNA methylation remodeling. This aspect of the circulating blood NK cell methylome before and after activation is in contrast to the significant IFNG demethylation required by naïve CD4 cells differentiating into IFNG competent Th1 polarized cells.45-47 The IFNG locus of naïve NK cells has thus been claimed to be constitutively demethylated.48 However, a recent detailed bisulfite sequencing study of NK cell IFNG loci indicated that a highly restricted region of IFNG becomes demethylated during cell activation and that this alteration confers a higher IFNG competence to fully differentiated NK cells.49 Our array experiments contained 2 CpG probes (cg00848007, cg26227465) located within this influential region in exon 1 and the immediate 5′ sequence. Targeted analysis of these specific probes confirmed significant demethylation at these loci following NK activation (data not shown), although the changes would not meet genome-wide significance. These results speak to the fact that subtle and highly localized epigenetic remodeling can be obscured when using only genome wide discovery.
Reprograming of DNA methylation within regulatory regions of transcription factor genes has been observed during immune activation. In murine CD8 cells following infection with the lymphocytic choriomeningitis virus (LCMV),50 methylation changes were evident at day 8 post-infection in cells demonstrating a fully differentiated effector CD8 cell phenotype. Our NK cells at day 9 post-activation share important methylation features with LCMV activated CD8 T cells including the demethylation of BHLHE40 (also known as Dec1, Stra13, Sharp2, or Bhlhb2). BHLHE40 belongs to a family of basic helix-loop-helix transcriptional regulators that is an attractive candidate for activation as it is induced in NK cells in response to cytokines (IL-2, IL-12, IL-15, INFG, TNFα) that enhance cytotoxic function.51 Additional roles of BHLHE40 and other differentially methylated transcription factors can be found in Table 3.13,52-57
Based on its role in T and NKT development and apparent significance in NK cell activation, several ddPCR assays for quantitating BHLHE40 demethylation were designed and evaluated. Cloning and bisulfite sequencing revealed broad and heterogeneous demethylation within the exon 5 region. Based on results, 2 CpG sites upstream of the array CpG probes were targeted for assay design. This upstream region contains a higher CpG density making it more favorable for MS-PCR design. This ddPCR assay accurately predicted the levels of unmethylated CpG sites in the exon 5 region. Consistent with its role in several lymphocyte lineages but not in myeloid cells, we observed low levels of BHLHE40 methylation in T, B and NK cells but high levels in granulocytes and monocytes. Among non-activated cells the highest prevalence of unmethylated BHLHE40 was found in naïve NK cells.
Among transcription factors, significant hypomethylation of the proximal ZBTB32 gene (also known as the Fanconi anemia zinc finger, or FAZF) and many fold induction of its expression stand out as notable features of NK cell activation in our studies. ZBTB32 is a member of the BTB-ZF transcription factor family, which can recruit chromatin-remodeling co-factors to regulate gene expression,58 and plays an essential cell-intrinsic and lineage-specific role in promoting NK cell expansion after viral infection.55 Recent studies focused on BTB-ZF family members as targets of epigenetic reprograming in CMV infection showed unmethylated CpGs in ZBTB32/FAPZ in HCMV memory NK cells but its expression was not elevated compared with naïve NK cells21 in contrast to our observed 600-fold increase in mRNA expression.
Among other functions, the NFAT family of transcription factors regulates cytokine gene expression by binding to the promoter/enhancer regions of antigen-responsive genes, usually in cooperation with heterologous DNA-binding partners. NFATC1 serves as regulator for many immune cell types such as T cells and natural killer cells.56 This study identified a locus in intron 9 that displayed significant, large demethylation (Δβ = −0.4; P =1.3 × 10−12; q = 5.6 × 10−8). The demethylation of intron 9 in activated NK deserves further study to understand what, if any, regulatory role it may have.
Activated NK cells contained a DMS in BACH2. BACH2 is a transcriptional repressor expressed by mature lymphoid cells that imposes regulation on humoral and cellular immunity. BACH2 helps maintain the balance between immune-stimulating and immune-regulating cells, or self-tolerance. Variants of BACH2 have been implicated in the susceptibility to several autoimmune diseases and allergy.59,60 Although a role for BACH2 in NK cells has not been determined, it plays an important role in lineage stabilization of Tregs.57 Given the significant commonalities between the cell types, BACH2 may function similarly in NK cells.
The CpG loci with the largest difference in DNA methylation upon NK activation occurred within the Mastermind-like 2 (MAML2) gene. MAML2 is a transcriptional co-activator of NOTCH,61 which is a highly conserved cell-cell communication pathway involved in numerous processes during embryogenesis, development, and hematopoiesis.62 Notch signaling enhances NK cell maturation and increases cytolytic effector capacity and cytokine secretion.63 Little information exists on the epigenetic features of MAML2. One study found the MAML2 promoter became hypermethylated in senescent mesenchymal stem cells after long term culture.64 We found another Notch related gene, the actin-binding protein involved in cytoskeletal organization,65 TAGLN3, demethylated in activated NK cells.
Significant hypomethylation of loci within the CISH gene was found in this study. CISH encodes a protein belonging to the cytokine-induced STAT inhibitor (CIS), also known as suppressor of cytokine signaling (SOCS) protein family. CIS family members are known to be cytokine-inducible negative regulators of cytokine signaling. Previous studies observed that DNA methylation levels within the second exon DMS of CISH was 100% in naïve CD4+ T cells, 52% in effector cells, and 13% in memory cells; an additional 17% of DMRs were hypermethylated in naïve cells, intermediately methylated in effector cells, and hypomethylated in memory cells.66
Since hypermethylation was so uncommonly observed in these activated NK cells, it is of interest to note RPS6KA1, ribosomal protein S6 kinase, polypeptide 1. Also known as p90RSK, this is a member of the ribosomal S6 protein serine/threonine kinases that regulate diverse cellular processes such as cell growth and differentiation. Putative methylation-based regulation of RPS6KA1 may result in modulation of the mTOR pathway, which has recently been recognized as important in the coordinate functioning of immune and metabolic activities in T cells.37
It is interesting to note that genomic regions previously identified in GWAS studies of immune related diseases co-localize with novel CpG loci differentially methylated in activated NK cells in the current study. The convergence of GWAS with EWAS of NK activation is an indirect but suggestive clue to the potential importance of DMRs in immune regulation.
A robust association of CLEC16A demethylation with NK activation points to another novel pathway previously unappreciated in NK biology. CpG site cg03995533 was significantly unmethylated in activated NK cells. This site at 16p13.13 lies in intron 18 of CLEC16A, a region identified in GWAS as a susceptibility locus for several autoimmune diseases including T1D and multiple sclerosis.67 CLEC16A function is largely unexplored but it is thought to play a role in regulating mitophagy, a selective form of autophagy necessary for mitochondrial quality control.68 CLEC16A is expressed in NK cells67 but is also abundant in CNS white matter and blood monocyte-derived dendritic cells from MS patients, in which it strongly co-localized with human leukocyte antigen class II. Alleles of the intron 18 GWAS SNP rs12708716, are associated with CLEC16A expression in pancreatic islet cells68 and thymic tissue. A very recent GWAS of common variable immunodeficiency disorder identified rs17806056, also in intron 18, as a significant susceptibility factor for that disease.69 This SNP lays 1,995 bp upstream of the activated NK DMS. The convergence of replicated GWAS signals with regions affected by differential DNA methylation during NK activation supports the plausible regulatory importance of these epigenetic changes.
Multiple CpGs within intron 1 and the proximal promoter regions of the NCF4 gene were hypomethylated in activated NK cells. The NCF4 protein is a cytosolic regulatory component of the superoxide-producing NADPH-oxidase whose function is most strongly associated with the antibactericidal action of reactive oxygen species produced by myeloid linage cells. In other immune cells, NADPH oxidase plays an important role at several levels in autophagy by regulating signaling, antigen processing and phagosomal pH.70,71 NCF4 deficiency leads to aberrant degradation and accumulation of the pro-inflammatory cytokine IL1B, which contributes to granulomatous colitis in inflammatory bowel disease.72 It is of interest that GWAS SNPs implicated in Crohn's disease include NCF4 with the most consistent associations (rs4821544) localizing to intron 1 of the gene.73 Our strong hypomethylated signal at cg02532700 of intron 1 is 849 bps downstream of rs4821544. Other unmethylated loci occurred nearby and in the proximal promoter region of the gene.
Tumor necrosis family α induced protein 3 (TNFAIP3) is a well-known negative regulator of the NF-kB activation pathway. TNFAIP3 can attenuate the NF-kB activity triggered by signaling from TNF and Toll-like receptors,74 and is implicated in several autoimmune diseases. In activated NK cells, we found the CpG 08919597 site consistently hypomethylated; this site appears to be highly conserved, located within a TF binding cluster and a very strong DNaseI site/enhancer in ENCODE experiments. It also lays approximately 1200 bp from GWAS SNP candidates for autoimmune disease. TNFAIP3 is considered a tumor suppressor inactivated in cancer through deletion, mutation or alternatively by DNA hypermethylation. TNFAIP3 hypermethylation occurs in MALT lymphoma but not in chronic lymphocytic leukemia nor in multiple myeloma.75 The tumor related DNA methylation occurs within 18 CpGs sites of the gene promoter in contrast to the intronic location that we identified.
Dramatic and stable epigenetic alterations within innate immune cells after infection by human pathogens have been documented. Long-lived “memory-like” NK cells have been identified that have distinct cell surface phenotype and functional properties. Recent studies implicate DNA methylation as an important mechanism leading to the cardinal features of these cells. Specifically, NK cells from HCMV-infected individuals include distinct subsets of memory-like NK cells that are deficient in multiple transcription factors (ZBTB16/PAZF) and signaling proteins, the latter including the tyrosine kinase SYK20 and the signaling adaptors FCERG1 and EAT2.21 These insights led us to ask if these CMV activation features were also present in our NKp46/CD2/IL2 activated NK cells. The analysis included 2 of the SYK CpG sites that were reported to be hypermethylated in memory NK cells following CMV infection.20 Β-values for these 2 probes were both only sparsely methylated; thus do not indicate significant differences in SYK methylation. This suggests hypermethylation at these sites may be unique to NK memory. Similarly FCERG1 and EAT2 loci implicated in HCMV were broadly hypomethylated and unaffected by NK activation. ZBTB16 was also hypomethylated and, in contrast, the ZBTB12 member was hypermethylated. These contrasting methylation features may help to distinguish NK cells activated through disparate exposures.
Our data also have important implications for the proper interpretation of epigenome-wide association studies carried out in whole blood. As demonstrated convincingly by Bauer et al. (2015),76 associations of small changes in DNA methylation of specific loci associated with environmental exposures can be explained by exposure-associated changes in the leukocyte composition. The authors show that smoking-associated alterations of a CD3 lymphocyte subtype explain the reports by multiple authors of tobacco-associated DNA methylation alterations at the GPR15 locus. Because it is possible that tobacco smoking activates NK cells, we examined whether the long list of loci reported to have altered methylation associated with smoking77 might also be found as DMSs indicative of NK activation. We observed that change in methylation of the AHRR locus is both a marker of NK activation and smoking78 and that several other DMSs associated with NK activation that we found here have been reported by others to be associated with smoking, including: MYO1G (cg07826859), NCF4 (cg02532700), and GPR68 (cg05875421). This adds evidence to support the need for cautious interpretation of EWAS data, suggesting that the underlying leukocyte biology is responsive to exposure and may represent the signal captured in the EWAS data.
Detailed understanding of the methylome of activated NK cells has implications for emerging therapeutics. NK cells are seen as promising agents for cell-based cancer therapies.79 However, there is a great need to develop more effective protocols to activate and expand high numbers of NK cells ex vivo to be able to infuse sufficient numbers of functional NK cells to cancer patients. The methylation profile established here provides an epigenetic template to evaluate such protocols. In addition, new classes of epigenetic drugs that target important transcriptional mechanisms are being developed that have high specificity and biological activity.80 For example, EP300, which was identified in differential methylated loci in both hyper- and hypo-methylated NK loci, is a bromodomain containing transactivator that has been successfully targeted with small molecule inhibitors.81 The methylation profile of activated NK cells described here would predict that bromodomain inhibitors targeting EP300 could modulate NK activation, which may be beneficial or undesirable depending on the specific application.
The Illumina 450K panel is extensive but only interrogates a fraction of existing CpG sites in the genome. The activation protocol we used is an artificial ex vivo approach that may not represent physiologically relevant in vivo immune reactions. Different activation protocols may lead to different modifications of the methylome, thus our profile may be one of many possible signatures. Also, the output of arrays may reflect the average of subsets within the cell population. Thus, the differential methylation signature of activated NK cells described here may be a composite reflecting the sum of several subsets. Our aim here was to search for potential markers that would faithfully reflect NK activation status (not predict protein levels); additional functional study of DMSs identified here will therefore be necessary to reveal the relationship between differential methylation and protein level and function.
Long-term culture and stimulation but not shorter-term activation times were associated with robust and reproducible changes in DNA methylation within human NK cells. At day 9 post-activation, most methylation changes consisted of hypomethylation that occurred predominately outside of CpG islands and instead featured gene bodies or CpG sparse regions. Genes previously identified as induced in activated lymphocytes, especially T cells, and that are regulated in part through DNA methylation, were prominently represented among the most significant differentially methylated sites in this study. Well-known immune activation loci included TNFA family members, CD96, Th2 cytokines (IL13) and co-stimulatory loci. We observed other novel loci including transcription factors previously identified in lymphocyte activation but not implicated as targets of DNA methylation (BHLHE40, ZBTB32, NFATC1). Genes implicated in specific haematopoietic motility (e.g., MYO1G) and Notch signaling (MAML2, TGLN3) were also identified for the first time. Also noteworthy were less well characterized genes implicated through genetic studies in autoimmune disease (CLEC16A, NCF4, BACH2) and for whom the relevant genetic susceptibility loci co-localized with the regions affected by DNA methylation changes in NK cells. Our ability to target BHLHE40 demethylation using ddPCR provides a potential new tool for NK studies. The detection of BHLHE40 demethylation in the whole blood of (healthy) individuals suggests this marker could be investigated in clinical populations. Potential applications include assessing immunomodulatory therapies or monitoring the levels of ex vivo activated NK cells in cell mediated therapies. Finally, the sum of many DMSs identified here could be used in population EWAS studies to identify and adjust for the contribution of activated NK cells as a strategy to avoid confounding by the effects of immune activation on DNA methylation profiles.
A discovery-replication experimental design was applied to guide our identification of valid CpG loci associated with NK cell activation (Suppl Fig. 1A). A strict false discovery threshold (q < 0.003) was applied to the 450K methylation data after filtering out those CpG loci that may have been affected by gender (X,Y chromosomes) or autosomal polymorphisms (SNP >0.05). The discovery set consisted of 6 different blood donors whose NK cells were studied at time 0 and day 9 post-activation, 4 of whom had additional samples taken at time 16 h. The replication phase consisted of a second set of 6 subjects processed and arrayed independently at a later time. Only loci that met the threshold criteria (q < 0.003) in both phases were included in the final data set.
Within 3–5 h of blood draw, freshly generated leukoreduction filters were used to isolate NK cells; the filters were supported on a clamp/ring-stand and back-flushed with ~250 ml of room temperature, sterile, PBS (without calcium or magnesium), pH 7.2, 2 mM EDTA (Miltenyi Rinsing Solution, cat. no. 130-091-222). Flushing was carried out with ˜50 ml at a rate of approximately 2 min per syringe of PBS into a single use sterile Erlenmeyer flask, while gently massaging the filter to enhance cell release. Flushed cell suspensions were subjected to Ficoll Paque gradients (GE cat. no. 17-1440-02, density 1.077 g/ml) in 50 ml conical tubes. The resulting buffy coat layers were collected and washed several times to remove platelets prior to the NK cell enrichment procedure.
Untouched NK cells were obtained using the Miltenyi NK Cell Isolation Kit (cat. no. 130-092-657) per manufacturer's instructions. Immediately following isolation, sub-samples were removed for both flow cytometry sample preparation and nucleic acids extraction. The remaining cells were activated using the Miltenyi NK Cell Activation & Expansion Kit (cat. no. 130-094-483), with beads prepared as directed by the manufacturer. This activation is a combination of CD2 and CD335 (NKp46) coated microbeads and exposure to 500 IU/ml IL2 (Miltenyi 130-097-744). The NK cells were cultured at a density of 1–1.5 × 106 cells/ml for 9 d under standard conditions (37°C, 5%CO2) in CellGenix SCGM (cat. no. 20802-0500) supplemented with 10% heat inactivated FBS (Life Technologies cat. no. 16000-044) and 1x L-glutamine/Pen/Strep mix (2 mM L-glut, 100 U/ml penicillin, 100 µg/ml streptomycin; Life Technologies cat. no. 10378016).
Fresh, magnetic-activated cell sorting (MACS) isolated and 9-day cultured NK cells were characterized by flow cytometric techniques. Samples were analyzed for viability, NK purity and activation status using a 5-stain panel (Suppl Fig. 2 and Suppl Table 1). Compensation controls were prepared using anti-mouse compensation beads (Life technologies ABC Anti-mouse Bead kit, cat. no. A10344). Heat killed cells (65°C, 1 min followed by 1 min on ice) were used as compensation controls for viability dye. All cells were washed twice and suspended in protein-free PBS. When IFN-γ was measured, an intracellular stain was employed. An aliquot of cells for IFN-γ staining were first incubated with a Golgi blocker, brefeldin A (eBioscience cat. no. 00-4506), for 4 h. Cells were stained for surface markers per manufacturer's instructions, and fixed in IC Fixation Buffer (eBioscience cat. no. 00-8222). Intracellular staining was carried out after permeabilizing cells (Permeabilization Buffer, eBioscience cat. no. 00-8333), as directed by the manufacturer. All fixed flow cytometry samples were stored at 4°C in the dark and analyzed within ~36 h on a FACS Aria (UCSF Laboratory for Cell Analysis) or FACS AriaUII (Brown University Flow Cytometry and Sorting Facility).
A killing assay was developed at Janssen Pharmaceutical and employed at Brown University to assess the killing phenotype of activated NK cells, using target Daudi cells. Killing was determined by measuring free GAPDH in the media/supernatant. Daudi cells (ATCC CCL-213) were cultured in RPMI 1640 supplemented with L-glutamine, 1% Pen-Strep and 10% heat inactivated FBS, and incubated at 37°C in a humidified incubator with 5% CO2. Cultures were maintained at a density between 3 × 105–3 × 106 cells/ml. Killing assay wells of U-bottom tissue culture plates were set up as follows: medium only (background control), Daudi cells only for target cell spontaneous lysis, NK cells only for effector cell spontaneous lysis, Daudi cells only for maximum target cell lysis via chemical disruption, effector:target (5:1), and effector:target (5:1) in the presence of 5 µg/ml Rituximab. Wells that included Rituximab were incubated with 40,000 viable Daudi cells in a final volume of 75 µL medium at room temperature for 15 min for opsonization. Plates were centrifuged (300 x g, 5 min) and medium removed. Washed, viable, day 9 activated NK cells (2 × 105) in complete medium were added to wells as described above. Cells were gently pelleted (300 x g, 1 min) to promote cell contact, and incubated at 37°C for 4 h for NK killing of Daudi cells. Following incubation, 10 µL of lysis buffer (Component 6, part #3035 from aCellaTox kit, Cell Technologies) was added to the subset of Daudi only wells designated for the maximum lysis control and incubated 15 min at room temperature. All wells were brought to a final volume of 200 µL with medium and centrifuge 750 x g, 1 min. Detection of GAPDH in the medium/ supernatant was carried out using reagents from the eBioscience GAPDH InstantOne ELISA kit (cat. no. 85-86131-11) with some procedural modification. Fifty microliters of experimental supernatant was added to 50 µL of 1:1 InstantOne detection and capture antibodies in the wells of the provided InstantOne ELISA assay plate. Plates were sealed and mixed (140–300 RPM) for 1 h at room temperature. Following incubation, the cell supernatant/antibody mixture was removed and plate wells washed 3 times with 200 µL /well wash buffer. One hundred microliters of detection reagent was added per well and incubated at room temperature for 10–15 min on an orbital shaker at 50 RPM. Reactions were stopped by the addition of 100 µL /well of kit stop reagent. Absorbance at 450 nm was read on a Molecular Devices Spectromax M2.
Specific NK killing percentages were calculated as follows:
Where experimental lysis is the signal from NK + Daudi cells, spontaneous lysis is the sum of signal from NK-only and Daudi-only wells, and maximum lysis is Daudi cells in the presence of chemical lysis solution.
NK cells were harvested at T = 0 (inactivated/naïve), T = 16 h (activated; discovery phase only), and T = 9 d (activated), and DNA extracted using the Norgen RNA/DNA Purification kit (cat. no. 48600) for discovery phase samples, and Zymo Universal Quick DNA Miniprep kit (cat. no. D4068) for replication phase samples. DNA was bisulfite treated and converted using the Zymo EZ DNA Methylation kit (cat. no. D5001) per manufacturer's instructions, including the alternative incubation conditions recommended for the Illumina InfiniumMethylation450 BeadChip array samples. Bisulfite converted DNA was sent to the UCSF Center for Human Genetics, Genomics Core facility for processing using Illumina chemistry, instruments and protocols.
Total RNA was purified from NK cell pellets using ALLPrep DNA/RNA Mini Kits (Qiagen cat. no. 80204) and cDNA created with random primers and SuperScript II Reverse Transcriptase (Invitrogen cat. no. 18604-022) based on vendor's instructions. Gene expression was assessed by comparative CT method (ΔΔCT) using IPO8, B2M, and GAPDH as the endogenous references genes. All expression assays were inventory Taqman Gene Expression Assays (LifeTechnologies; Suppl Table 2). Assays were carried out per manufacturer's instructions, using 2 µl of cDNA per 20 µl PCR. Reactions were run in triplicate on an Applied Biosystems 7300. ΔΔCT was calculated using 9 d activated samples as the “test” and T = 0 (freshly isolated NK cells) as the “calibrator.”
Genomic DNA for bisulfite sequencing was extracted from NK cells using ALLPrep DNA/RNA Mini Kits (Qiagen cat. no. 80204) and bisulfite treated with EZ DNA Methylation Kit (Zymo cat. no. D5001).
Bisulfite converted DNA was amplified with bisulfite-specific and methylation insensitive primers. PCR products were purified using either a Gel Extraction Kit (Qiagen cat. no. 28604) or PCR Purification Kit (Qiagen cat no 28104), and PCR products ligated into plasmid vectors using the TOPO TA Cloning system (Invitrogen cat. no. K4500-0). Plasmid-transformed Escherichia coli were cultured overnight on Luria-Bertani agar plates and the positive colonies were cultured overnight in Luria-Bertani media. Plasmid DNA containing PCR products were isolated with QuickLyse MiniPrep Kits (Qiagen cat. no. 27405) and 20 colonies from each experiment group were sequenced in both directions using an ABI 373 automated sequencer with dye primer chemistry and standard M13 primers.
As described in Suppl Fig. 1B, the most significant loci from the analysis were further prioritized for PCR assay development and ranked using 3 additional features: number of statistically significant probes within a gene, effect size and CpG density. Gene loci were scored based on the presence of multiple probe hits within a gene and the effect size (>0.3 β-value difference between activated and non-activated conditions), as well as the density of CpGs surrounding the top hit probe. Based on these feasibility criteria, several candidate genes were selected and tested for potential quantitative methylation assay development (data not shown). BHLHE40 was chosen for further exploration as a putative biomarker of NK activation.
Four methylation-specific PCR (MSP) designs were selected for BHLHE40 and tested using the QX100 Droplet Digital PCR system (Bio-Rad)(data not shown). The following assay was selected: forward primer 5′- TGGAGGTGGGTGATAAGTTG-3′, reverse primer 5′- ACTCAAATTTCCAAACATATACCC-3′, and probe 6FAM-5′- AAACACAAAAACACTCAAAACCTAA-3′-NFQ. Master mixes for design testing included: Droplet PCR Supermix (Bio-Rad, cat. no. 186-3024), primers and probe (Applied Biosystems), and EB buffer (Qiagen cat. no. 19086). The so-called “C-Less” reaction (forward primer 5′-TTGTATGTATGTGAGTGTGGGAGAGA-3′, reverse primer 5′- TTTCTTCCACCCCTTCTCTTCC-3′ and probe 5′- 6FAM-CTCCCCCTCTAACTCTAT-NFQ-3)’ was used to quantify the total input of bisulfite converted DNA into each sample and BHLHE40 copy number expressed as a ratio to the total number of copies.82
BHLHE40 and C-Less total input master mixes were aliquoted into 96-well plates and PCR reactions performed in monoplex using 80 ng of bisulfite converted DNA for all reactions. Samples were transferred one column at a time to an 8 channel Droplet Generator Cartridge (Bio-Rad, cat no 186-3008). QX100 Droplet Generator Oil (Bio-Rad, cat no 186-3005) was added, and droplets generated using the QX100 Droplet Generator (Bio-Rad). Typically 10,000 to 14,000 droplets were generated for each sample and transferred to an Eppendorf 96-well PCR plate (cat no 951020362), foil heat sealed, and placed into a S1000 thermal cycler (Bio-Rad). The PCR cycling conditions for gradient temperature testing were as follows: 95°C for 10 min, 40 cycles of 94°C for 30 sec and 63°C-56°C for 1 min, and ending with 98°C for 10 min and 4°C hold. Following PCR, the plate was transferred to the QX100 Droplet Reader (Bio-Rad) for droplet analysis/determination of either positive or negative for fluorescence/amplification. Data was analyzed using the QuantaSoft software from BioRad.
BHLHE40 ddPCR assay performance was examined using activated NK cells as a positive standard of demethylated exon 5, which was cloned and bisulfite sequenced. Samples of activated NK DNA were diluted with hypermethylated DNA from granulocytes to simulate the detection of BHLHE40 within a background of excess hypermethylated DNA. Dilution of activated NK DNA with the background copies observed in completely hypermethylated DNA was used to estimate the assay limit of quantification and dynamic range (data not shown).
Gene-associated hits from the overlapping loci (discovery and replication data sets) were analyzed for functional enrichment as well as pathway association. Functional enrichment was evaluated using the National Institute of Allergy and Infectious Disease DAVID Bioinformatics Database v6.7 http://david.abcc.ncifcrf.gov/.83,84 P-value (or EASE score) was set at 0.05 (please refer to https://david.ncifcrf.gov/helps/functional_annotation.html#E3 for more detailed explanation). The lower the P-value (or EASE score), the more enriched the gene terms in the dataset compared with the representation of that gene-term in the human genome. Reference data was limited to using only human genes found in the selected (“Official Gene Symbol” setting in DAVID) annotation databases. Pathway analyses (“canonical pathway” module) were run on Ingenuity Pathway Analysis software (Qiagen). Analyses were run using a stringent filter for “human,” which returned results most relevant to humans. Additionally, analyses were run using all tissue and cell types, then repeated constraining the input/reference to immune cell information only. Non-coding associated loci were evaluated in a similar manner using the open access program GREAT (version 3.0.0; Suppl table 7). This program was designed to interrogate numerous annotation/ontology databases to calculate enrichment for proximal and distal, non-coding genomic regions. Analysis was carried out using the program default association rule parameters.
The significant differentially methylated loci (385) were also assessed in the context of transcription factor binding site. The TFBS for 161 transcription factors from the ENCODE database was interrogated using the “Txn Factor ChIP” track and wgEncodeRegTfbsClusteredV3 table through the UCSC Genome Browser (hg19 assembly) site.
Methylation signals from Illumina 450K methylation arrays were preprocessed using the functional normalization method “preprocessFunnorm” in the minfi package85 86 to control for unwanted variation among plates and samples. To avoid interference from polymorphisms affecting probe hybridization, all probes containing single nucleotide polymorphisms (SNPs) were removed as were probes contained on the X- and Y-chromosomes. Additionally, cross-reactive probes were removed using criteria from Chen, et al. 2013.87 Models for time 0 vs. 16 h activation, and 0 vs. 9 d activation found F statistics and P-values using an empirical Bayes method. These calculations were made using the function “eBayes” in the “limma” package.88 False discovery rate (FDR) was corrected for using the Benjamini-Hochberg adjustment for multiple testing.89 CpG loci were ordered by FDR-adjusted q-value, and the top 1000 loci selected for further consideration. Subsequent analyses focused on gene-associated CpG loci.
No potential conflicts of interest were disclosed.
We acknowledge the helpful suggestions and discussions of Dr. Shichun Zheng.
JKW was supported by the Robert Magnin Newman endowment in Neuro-oncology. This work was supported by Janssen Pharmaceuticals.