Selective DNA Demethylation Accompanies T Cell Homeostatic Proliferation and Gene Regulation in Lupus-Prone lpr Mice

Systemic lupus erythematosus (SLE) is characterized by increased DNA demethylation in T cells, although it is unclear whether this occurs primarily in a subset of SLE T cells. The process driving the DNA demethylation and the consequences on overall gene expression are also poorly understood and whether this represents a secondary consequence of SLE or a primary contributing factor. Lupus-prone lpr mice accumulate large numbers of T cells with age because of a mutation in Fas (CD95). The accumulating T cells include an unusual population of CD4−CD8−TCR-αβ+ (DN) T cells that arise from CD8+ precursors and are also found in human SLE. We have previously observed that T cell accumulation in lpr mice is due to dysregulation of T cell homeostatic proliferation, which parallels an increased expression of numerous genes in the DN subset, including several proinflammatory molecules and checkpoint blockers. We thus determined the DNA methylome in lpr DN T cells compared with their CD8+ precursors. Our findings show that DN T cells manifest discrete sites of extensive demethylation throughout the genome, and these sites correspond to the location of a large proportion of the upregulated genes. Thus, dysregulated homeostatic proliferation in lpr mice and consequent epigenetic alterations may be a contributing factor to lupus pathogenesis.


INTRODUCTION
T cells from patients with systemic lupus erythematosus (SLE) are known to manifest evidence of activation and autoreactivity (1,2). They also contain increased levels of DNA demethylation, one of the main epigenetic regulators of gene expression (3)(4)(5). In addition, certain medications that are known to promote DNA demethylation, such as hydralazine and procainamide, can provoke autoreactivity of T cells and drug-induced lupus (6). However, the mechanism driving DNA demethylation in SLE remains obscure.
This results in the accumulation of T cells that would ordinarily undergo programmed cell death during homeostatic proliferation (8). Among the accumulating T cells is a subset of polyclonal CD4 − CD8 − TCR-αβ + (DN) T cells that derive from CD8 + precursors during homeostatic proliferation (9,10). This subset is also present in human SLE and derives from CD8 + T cells (11)(12)(13). We have previously observed that compared with their CD8 + precursors, DN T cells from lpr mice have upregulated gene expression of numerous immune modulating molecules, including the cytolytic machinery of Fas-ligand, Granzyme B, and perforin, as well as inhibitory molecules such as PD-1, Lag3, and IL-10 (14). Initial analysis of one of these genes, Pdcd1 (PD-1), which is known to be regulated by DNA methylation (15), revealed extensive demethylation of the 5ʹ regulatory region in DN T cells compared with the CD8 + precursors subset (14). Based on these observations, we considered that DNA demethylation may occur more extensively in the genome of lpr DN T cells as part of the process of homeostatic proliferation. This might serve in part to explain the particular constellation of genes upregulated in these cells.

Mice
Mice were bred and housed in the Association for Assessment and Accreditation of Laboratory Animal Care International-approved animal facilities of The University of Vermont. Original breeding pairs of B6.MRL-Fas lpr /J (Fas lpr/lpr ) mice were obtained from The Jackson Laboratory (Bar Harbor, ME). All mice in these studies were on a C57BL/6 background and were used between 10 and 13 wk of age. All animal studies were conducted in accordance with the policies of The University of Vermont's Animal Care and Use Committee.

Reduced representation bisulfite sequencing
Genomic DNA was prepared using the Quick gDNA Micro Prep Kit (Zymo Research) according to the manufacturer's instructions. Bisulfite sequencing libraries were generated as described previously using a custom adapter-primer combination (16,17). Briefly, genomic DNA was digested separately with MspI and TaqI to enrich for CpG containing DNA and combined with methylated PhiX control DNA (New England Biolabs). The resulting DNA was used as input for the KAPA HyperPrep Kit (Roche Diagnostics) with the following modifications. Adapters for ligation contained fully methylated CpGs. Following ligation, adapter-ligated DNA was bisulfite treated using the EpiTect Bisulfite Kit (QIAGEN), and library amplification was performed using the KAPA HiFi HotStart Uracil+ polymerase (Roche Diagnostics) with custom indexing primers. Final libraries were quality checked on a Bioanalyzer 2100 (Agilent Technologies), pooled at equimolar ratio, and sequenced on an Illumina HiSeq1500/2500 RapidRun using a paired-end 2 × 75 Rapid Run flow cell at the Vermont Integrated Genomic Resource Core.

Reduced representation bisulfite sequencing analysis
The FASTQ files were quality checked using the Fastx Toolkit v0.0.14, adapter content trimmed using Cutadapt v1.12 (18), and data mapped to the mm10 version of the mouse genome using Bismark v0.13 (19). CpG methylation calls were computed using a custom R pipeline that is available upon request. Bisulfite conversion was assessed using the PhiXmethylated spike in DNA, which showed >99% conversion efficiency for all samples. CpG covered with at least 10 reads in each group were used for all downstream analysis. Differentially methylated loci (DML) were determined using the DSS package (20) and CpG with a false discovery rate-corrected p value <0.05 and <20% change in methylation between groups were considered significant. CpG were annotated to the nearest gene transcription start site. To identify transcription factor binding motifs at demethylated DML the HOMER findMotifs-Genome.pl script was used with the -size 400 setting. All data visualization and downstream analysis were performed using custom R scripts that are available upon request. DNA methylation data are available from the Gene Expression Omnibus under accession https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi? acc=GSE155293.

Meta-analysis with microarray data
To determine to what degree the most demethylated genes were part of the set of genes upregulated in DN T cells compared with CD8 + T cells, DNA methylation data were integrated with our previous microarray analysis (14). A total of 968 genes were both upregulated and demethylated DN T cells. With this set, Gene Ontology (GO) term enrichment and pathway analysis was conducted using Partek Genomics Suite, version 7.18 (Partek, St. Louis, MO). In addition, functional cluster and pathway analysis was performed using both National Institutes of Health (NIH) Database for Annotation, Visualization, and Integrated Discovery (https://david.ncifcrf.gov) and Ingenuity Pathways Analysis (Ingenuity Systems, www.ingenuity.com). The DML that mapped to the 968 upregulated genes was selected, resulting in 7649 CpG. Among those, 7267 DML were demethylated in the DN subset compared with the CD8 + subset.

Meta-analysis with assay for transposase-accessible chromatin sequencing data
To compute the overlap of DML and accessible chromatin regions, we analyzed previously defined assay for transposase-accessible chromatin sequencing (ATAC-seq) data (accession no. GSE83081) from unstimulated splenic naive CD8 + T cells and splenic effector CD8 + T cells at day 8 following lymphocytic choriomeningitis virus Armstrong infection (day 8 effector) that were positive for the gp33, gp276, and np396 tetramers (21). A set of accessible regions was assembled by merging the peaks from naive and day 8 effector samples using the HOMER mergePeaks function and converted to a bed file using pos2bed.pl (22). Because the ATAC-seq data were mapped to the mm9 mouse genome, all DML coordinates were converted mm9 using the UCSC Genome Browser liftover tool (23). The overlap of DML and ATAC-seq regions was computed using the bedtools window function with -w as the indicated distance (24). To determine overlap significance, the ATAC-seq regions were randomly shuffled across the genome using the bedtools shuffle command 1000 times and overlap with DML calculated as above for each distance window.
The p values were calculated as the number of times the permuted overlap was greater than the observed overlap divided by the number of permutations (1000), with p values <0.001 resulting from zero permutations being greater than the observed overlap. For plotting ATAC-seq data, the mm9 bigWig tracks were converted to mm10 using the bwtool (25) and the mm9ToMm10.over.chain.gz UCSC Genome Browser chain file.

Discrete genome demethylation of lpr DN T cells
DNA from DN T cells and the precursor CD8 + T cells of B6-lpr mice were subjected to reduced representation bisulfite sequencing (RRBS) (17). B6-lpr mice were chosen over MRL-lpr mice, as B6-lpr mice develop only minimal if any lupus manifestations, allowing better separation of the epigenetic findings from potential confounding processes driven by autoimmune disease. Three separate purifications of both Tcell subsets were made using a pool of three mice in each purification. A total of 2,350,129 CpG with a 10× coverage across both groups of T cells was used for downstream analysis, representing 11% of the roughly 22 million murine CpG. Analysis of the mean DNA methylation levels in each subset revealed that the DN subset had ~2.5% greater demethylation globally compared with the CD8 + T cell precursors (Fig. 1A, 1B).
Analysis of significant DML revealed 56,903 CpG (2.4% of total covered) that displayed >20% change in DNA methylation with a false discovery rate <0.05 (Supplemental Table I). Intriguingly, 54,500 (96%) of the DML between these two T cell subsets reflected a loss of DNA methylation in the DN subset, which was highly consistent and significant across the three separate experiments (Fig. 1C, 1D). Consistent with the change in distribution, many of the CpG were close to 100% methylated in all CD8 + T cell alleles and shifted to lowintermediate levels of methylation in DN T cells. Of particular interest was that multiple DML clustered at discrete genomic regions, suggesting the changes in DNA methylation were not random. For example, significant changes in DNA methylation along chromosome 1 cluster when viewed at different base pair resolutions (Fig. 1E). Findings for all chromosomes are provided in Gene Expression Omnibus under accession GSE155293. These data demonstrate that DN T cells display wide-ranging DNA demethylation but at discrete loci.

Overlap of demethylation sites and upregulated genes in DN T cells
We previously performed gene expression profiling for DN T cells and CD8 + T cells from B6-lpr mice and noted 1646 genes that were significantly upregulated in DN T cells (14). Among the upregulated genes was Pdcd1, and the finding that this locus was also demethylated suggested there might be a wider correlation of DNA methylation changes and gene expression of other upregulated genes in DN T cells. Therefore, DML were assigned to genes by annotating to the nearest transcriptional start site. Strikingly, 968 of the 1646 upregulated genes mapped to DML that lost DNA methylation in the DN subset ( Fig. 2A). Given that many of these genes are upregulated in DN T cells, the change in gene expression was correlated with the change in DNA methylation for each DML that mapped to a gene with significant gene expression changes. Consistent with a repressive role for DNA methylation, this analysis revealed an inverse relationship between the two datasets, with the majority of the demethylated loci mapping to genes that gained gene expression (Fig. 2B, top left quadrant).
A GO term enrichment and Kyoto Encyclopedia of Genes and Genomes pathway analysis was performed using this set of 968 upregulated genes to determine if there was an enrichment for specific cellular functions. Among the enriched GO terms were molecular mechanisms of cancer, PI3K signaling, PTEN signaling, and cell death and survival (Table  I). The leading Kyoto Encyclopedia of Genes and Genomes pathway was TCR signaling (Supplemental Table I). This was of interest given that the DN T cells arise from CD8 + precursors during repeated rounds of homeostatic proliferation, which requires recurrent TCR stimulation by autologous MHC/peptide complexes (26). Additionally, the upregulated genes were ranked according to the number of demethylated DML mapping to each gene (Supplemental Table II). The gene with the most DML and highest ranking was Foxp1, which is a known regulator of quiescence in T cells (27,28). This is consistent with earlier observations that lpr DN T cells are small senescent cells that do not proliferate when activated in vitro (29,30).

DML occur near regions of accessible chromatin in effector CD8 + T cells
Analysis of the gene expression changes for DNA methyltransferases and demethylases revealed that only Dnmt3b was significantly different (Table II), suggesting other mechanisms might explain the observed DNA methylation changes. Analogous to the development of DN T cells from CD8 + precursors during homeostatic proliferation, the differentiation of CD8 + T cells from naive to effector cells involves extensive remodeling of accessible chromatin and DNA methylation at cis-regulatory regions (21,31). The observed clustering of DML and their proximity to genes that were remodeled during effector CD8 + T cell differentiation suggested that the DNA methylation changes may be associated with similar cis-regulatory elements. Therefore, using a range of distance windows, the overlap of each DML to a region of accessible chromatin in naive or day 8 effector CD8 + T cells responding to lymphocytic choriomeningitis virus from a previous study (21) was calculated. An increasing overlap of DML and accessible regions was observed, which was greater than a set of randomly shuffled sequences, with 16% within 500 bp and 54% within 10 kb (Fig. 4A). Analysis of the accessibility in the 4 kb surrounding DML demonstrated a significant gain in accessibility in day 8 effector CD8 + T cells over naive CD8 + T cells (Fig.  4B, 4C). For example, the Il10 and Pdcd1 (PD-1) loci contained DML that both overlapped and occurred in proximity to regions that gain accessibility in day 8 effector CD8 + T cells (Fig. 4D). Thus, DNA methylation occurs in proximity to cis-regulatory elements that gain accessibility in effector CD8 + T cells.

DML are enriched for AP-1, T-BET, and EGR transcription factor binding motifs
Regions that demonstrate dynamic DNA methylation during CD8 + T cell differentiation are enriched for transcription factor binding motifs that play important roles in CD8 + T cell fate and function (31). To determine if similar mechanisms occurred at the demethylated CpG in DN T cells, the surrounding 200 bp of each DML was searched for enriched transcription factor binding motifs using HOMER (22). As expected, the top two scoring enriched motifs contained sites for the MspI (CCGG) and TaqI (TCGA) restriction enzymes used in the RRBS assay (Fig. 5A). The next top scoring motifs were for the transcription factors AP-1, T-BET, and EGR2, all of which are involved in CD8 + T cell effector function (32)(33)(34). To further examine the epigenetic changes, the levels of DNA methylation and accessibility were computed for the 400 bp surrounding both AP-1 and T-BET motifs. Both motifs demonstrated a loss in DNA methylation in DN T cells versus CD8 + T cells (Fig. 5B). Consistent with this finding, each motif showed a higher level of accessibility in day 8 effector compared with naive CD8 + T cells. These data suggest that DML in proximity to transcription factor binding sites drive an effector program as CD8 + T cells transition to DN T cells during homeostatic proliferation.

DISCUSSION
The current findings reveal an epigenetic program of global, yet highly selective, DNA demethylation accompanied by upregulation of numerous genes, both of which accompany T cell homeostatic proliferation in lpr mice. Unlike MRL-lpr mice, B6-lpr mice develop very little autoimmune disease with age. We thus intentionally performed these studies using B6lpr mice to separate the epigenetic changes that accompany T cell homeostatic proliferation from confounding factors that might be secondary to autoimmune disease. The findings therefore may be more broadly applicable to T cell homeostatic proliferation also in wildtype mice. In fact, we have previously observed that successive rounds of T cell homeostatic proliferation in B6 wild-type mice parallel gene expression changes nearly identical to B6lpr T cells, and a DN subset is also observed in B6 wild-type mice (14). Conceivably, epigenetic regulation of a similar constellation of genes may occur in wild-type mice during T cell homeostatic proliferation but are enhanced in lupus-prone mice.
T cells from SLE patients bear demethylated DNA (3)(4)(5), but the underlying cellular process regulating this was not fully defined. It was also not clear whether the demethylation occurred primarily in a subset of SLE T cells or globally in all T cells. Most of the studies on DNA methylation in lupus T cells have been conducted on CD4 + T cells, whereas the current study examined DNT cells and their precursor CD8 + Tcells. Nonetheless, there are remarkable similarities in upregulated gene expression between SLE CD4 + T cells and lpr DN T cells, including ITGAL (CD11a), CD40LG (CD40-ligand), IFNG (IFN-γ), RFR1 (perforin), GADD45A, and CXCR3 (5), suggesting a possible common mechanism.
Because DN T cells also occur in human SLE (2,11), conceivably this subset in human SLE may also bear a substantial portion of the DNA demethylation in SLE T cells. Although most humans with SLE do not bear Fas gene mutations, it is possible that mechanisms operative in SLE might accelerate T cell homeostatic proliferation, resulting in a similar epigenetic program.
Recently, a study of PBMC from 17 monozygotic and dizygotic twin pairs discordant for SLE revealed extensive demethylation at 807 CpG sites corresponding to 49 genes in the affected twin compared with their healthy twin (35). This was not observed in twinsdiscordant forrheumatoid arthritis or diabetes mellitus. This study also found that the SLE patients had reduced mRNAlevels of the DNA methyltransferase, Dnmt3b, which we also observed to be decreased in lpr DN T cells (Table II). A further similarity is that several of the immune genes that were both demethylated and upregulated in the SLE twin study were also demethylated and upregulated in lpr DN T cells, including Il10, Grb10, Gfi1, Padi4, Cd9, and Aim2 (14,35). In addition, DNA demethylation of a particular gene, Tnfsf7 (CD70), has been observed in 16-wk autoimmune MRL-lpr mice compared with 5-wk mice prior to onset of autoimmunity (36). This was associated with a reduction in levels of DNMT1 (36).
The current findings also have striking similarities to our recent epigenetic and transcriptional analyses of B cells in human SLE (16). A subpopulation of IgD − CD27 − CXCR7 − CD11c + (DN2) B cells is expanded in SLE and has been linked to disease (37). DN2 B cells also share some similarity with autoantibody-associated B cells described in mice (18,38). Interestingly, the AP-1, T-BET, and EGR transcription factor binding motifs that were observed in this study are also enriched in the SLE-specific accessible chromatin (16). Lpr DN T cells express high levels of the AP-1 complex (39) and manifest high levels of Ifng (14). Additionally, the gene expression profiling showed a parallel dysregulation of similar gene networks, among them regulation of cell cycling, glycolysis, oxidative phosphorylation, and apoptosis. This parallels the known high levels of both glycolysis and oxygen consumption by lpr DN T cells (40,41). Additional common groups of upregulated genes included TCR signaling and protein phosphorylation. This is also consistent with the known constitutive phosphorylation of many signaling proteins in lpr DN Tcells, including several TCR signaling components, such as CD3ζ and Fyn (42,43). Scharer

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript In addition to the numerous upregulated genes related to cell cycling and cytokine production, both lpr DN T cells and SLE DN2 B cells also manifest upregulation of the cell cycle checkpoint blocker PD-1 (Pdcd1) and IL-10, and the sites of both genes correspond to sites of demethylation (16). It is of interest that homeostatically repopulating B cells following depletion with rituximab also express high levels of IL-10 (44). Collectively, these findings may suggest that lpr DN T cells, SLE DN2 B cells, and B cells following rituximab share a common epigenetic and transcriptional program linked to increased homeostatic proliferation.
A recent study showed that treatment of MRL-lpr mice with 5-azacytidine (5-Aza), a chemical analogue of cytidine that inhibits DNA methylation, enhanced autoimmunity, whereas targeting 5-Aza to only CD4+ or CD8 T cells alleviated disease (45). The reason for this seeming paradox was not fully explored, but 5-Aza also inhibits DNA replication, which could reduce the adenopathy of these mice when targeted to T cells, a feature not described in the study.
There may be several reasons for the progressive demethylation in lpr DN T cells.
Expression of the DNA methyltransferase, Dnmt3b, is reduced in this subset compared with the CD8 + precursors, which could contribute to this genotype. An alternative is that rapidly proliferating cells can exceed the ability of DNA methyltransferases to keep pace with the rapid DNA replication rate (46). Consistent with this, we have previously observed that the lpr DN T cells have undergone very rapid proliferation in vivo, with up to 18% replicating during a single 24-h period, as defined by BrdU uptake (14,26). An additional possibility is oxidative stress-induced demethylation. T cells from patients with active lupus manifest increased mitochondrial oxidative phosphorylation and reactive oxygen species (ROS) (47,48). Oxidative stress of CD4 + T cells results in decreased levels of DNMT1, DNA demethylation, and upregulation in expression of several genes (4). Moreover, adoptive transfer of oxidant-treated CD4 + T cells into syngeneic mice caused anti-dsDNA Ab and glomerulonephritis (49). Consistent with this, as noted earlier, lpr DN T cells manifest increased oxygen consumption compared with CD8 + precursor cells in part because of increased mitochondrial mass that parallels homeostatic proliferation (41).
Several of the demethylated and upregulated genes in lpr DN T cells are associated with inflammation (FasL, GzmB, Prf1, and Ifng) and immune exhaustion (Pdcd1 and Lag3). This might help explain the clinical immunology paradox of individuals with immunodeficiency syndromes either genetically, from chemotherapy, or because of HIV, which suddenly develop autoimmune syndromes. A dramatic example is the sudden onset of psoriasis and psoriatic arthritis in HIV + individuals (50). Conceivably, the lymphopenia in these conditions could lead to accelerated T cell homeostatic proliferation with resulting upregulation of inflammatory molecules. Conversely, SLE patients, bearing a seemingly overactive immune system, are nonetheless prone to infections and often respond poorly to vaccinations (51). It is possible that upregulation of PD-1 and Lag3 in SLE T cells renders them less responsive to new activation.
In summary, homeostatic proliferation of T cells manifests a broad program of both genetic and metabolic changes that could influence immune function and inflammatory autoimmune conditions. This includes the increased mitochondrial mass and size in T cells undergoing homeostatic proliferation (41). This contributes to high oxygen consumption rates, ROS production that can induce oligomerization of MAVS, and increased type I IFN (41,52). ROS and cell proliferation may also contribute to DNA demethylation (4). The current studies, thus, expand our knowledge of functional modifications during T cell homeostatic proliferation to reveal an epigenetic program of DNA demethylation at selective sites throughout the genome, contributing to upregulation of several immune response genes.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.