Epigenetic alteration contributes to the transcriptional reprogramming in T-cell prolymphocytic leukemia

T cell prolymphocytic leukemia (T-PLL) is a rare disease with aggressive clinical course. Cytogenetic analysis, whole-exome and whole-genome sequencing have identified primary structural alterations in T-PLL, including inversion, translocation and copy number variation. Recurrent somatic mutations were also identified in genes encoding chromatin regulators and those in the JAK-STAT signaling pathway. Epigenetic alterations are the hallmark of many cancers. However, genome-wide epigenomic profiles have not been reported in T-PLL, limiting the mechanistic study of its carcinogenesis. We hypothesize epigenetic mechanisms also play a key role in T-PLL pathogenesis. To systematically test this hypothesis, we generated genome-wide maps of regulatory regions using H3K4me3 and H3K27ac ChIP-seq, as well as RNA-seq data in both T-PLL patients and healthy individuals. We found that genes down-regulated in T-PLL are mainly associated with defense response, immune system or adaptive immune response, while up-regulated genes are enriched in developmental process, as well as WNT signaling pathway with crucial roles in cell fate decision. In particular, our analysis revealed a global alteration of regulatory landscape in T-PLL, with differential peaks highly enriched for binding motifs of immune related transcription factors, supporting the epigenetic regulation of oncogenes and genes involved in DNA damage response and T-cell activation. Together, our work reveals a causal role of epigenetic dysregulation in T-PLL.


T-PLL patients and healthy subjects. The Mayo Clinic patients who had a clinical diagnosis of T-PLL
and donated their blood samples to Mayo Clinic chronic lymphocytic leukemia (Mayo Clinic IRB 1827-00) and lymphoma (Mayo Clinic IRB 118) tissue banks were used in this study. After screening all T-PLL samples in the above tissue banks, only six were stored with more than 50 million fresh peripheral blood mononuclear cells (PBMCs) that were used for ChIP-seq in this study ( Table 1). Four of the six samples were used for total RNA sequencing; two of the six samples plus another three samples were used for mRNA sequencing. All but one of the nine T-PLL cases had the characteristic rearrangement involving TCL1A. The Mayo Clinic blood bank provided peripheral blood samples from three age-matched healthy individuals for ChIP-seq and total RNA sequencing and another four for mRNA sequencing (Mayo Clinic IRB 07-05543). Samples were collected with written consent and approval from the institutional review board at Mayo Clinic. We used PBMCs from the Table 1. Clinical information about the nine T-PLL cases. *Treatment history prior to sample collection. Patient P1 received a three-month romidepsin treatment, which was discontinued one month prior to sample collection. + Time from diagnosis to therapy. 46,X,-X,t(6;6)(q27;q15),inv (14) (q13q24), + mar [15]/ RNA-seq data analysis. Raw reads were aligned to the hg19 reference genome using the spliced transcripts alignment to a reference program (STAR, v2.5.2b, https:// github. com/ alexd obin/ STAR) 8 Read counts per gene were estimated for the Ensembl gene annotation, using the featureCounts tool in the Subread package (v1.5.1, https:// subre ad. sourc eforge. net) 9 . Gene expression was quantified as reads per kilobase of exon per million mapped reads (RPKM). Protein-coding genes differentially expressed between T-PLL and normal were identified using the edgeR package (v3.16.5, https:// bioco nduct or. org/ packa ges/ edgeR) 10 at the cutoff of 5% false discovery rate (FDR) and threefold change. Briefly, genes without reads in all libraries were excluded; for the remaining genes, raw read counts were normalized using the trimmed mean of M-values (TMM) method within edgeR. We used conditional maximum likelihood method within edgeR to estimate common negative binomial dispersion across genes, which is proved to be effective for small sample size.
To functionally categorize the differentially expressed genes, pathway and gene set enrichment analyses were done separately for the genes up-and down-regulated in T-PLL. For pathway analysis, we used the ReactomePA package (v1.18.1, https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ React omePA. html) 11 that takes the Reactome database 12 . For gene set enrichment analysis, we used the GOstats package (v2.40.0, https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ GOsta ts. html) 13 , which calculates p value for the overrepresentation of gene sets in Gene Ontology (GO) Biological Process and Molecular Function.
To highlight the gene expression changes between T-PLL and normal and the level of heterogeneity within T-PLL, we performed unsupervised clustering of expression levels for the top 10,000 protein-coding genes, using the pheatmap package (v1.0.8, https:// cran.r-proje ct. org/ web/ packa ges/ pheat map/) in R. Gene expression RPKM was log 2 transformed and quantile normalization using limma package (v 3.30.13, https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ limma. html). The 10,000 genes were selected from autosomes that had RPKM values ≥ 0.5 in at least one of the samples and showed the largest variation across samples. Total RNA and mRNA sequencing data were separated when used in the above differential analysis and unsupervised clustering.
ChIP-seq data analysis. ChIP-seq data were analyzed with the HiChIP pipeline with minor modifications 14 .
To perform unsupervised clustering for H3K4me3 and H3K27ac peaks, we first merged peaks if they were overlapped by at least 1 bp. For each merged peak in each sample, we estimated the number of mapped pairs in the IP and input using BEDTools. Read counts in IP were input-subtracted, normalized to 10 M mapped reads, log 2 transformed and quantile normalized. We retained merged peaks that were from autosomal chromosomes and present in at least two samples. For H3K27ac, we focused on merged peaks that are located in the enhancers and removed those that overlapped the TSS ± 2.5 kb regions from protein-coding genes, as described 20,21 . The top 20,000 merged peaks with the largest variation across samples were retained for unsupervised clustering. H3K4me3 peaks were processed similarly, except that we focused on those in the promoters of proteincoding genes, defined as the TSS ± 2 kb regions 22  www.nature.com/scientificreports/ down-regulated in T-PLL, we generated aggregate per-million signal density profiles over TSS ± 2 kb using ngs. plot tool (v2.63, https:// github. com/ shenl ab-sinai/ ngspl ot) 23 . We used the DiffBind R package (v2.14.0, https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ DiffB ind. html) 24 to identify the changes of H3K27ac and H3K4me3 between the six T-PLL and three normal subjects, and between the four T-PLL samples and three normal controls that also have total RNA sequencing data. Only merged peaks that were present in ≥ 2 samples were included. Data were TMM normalized and differential peaks were identified using the edgeR method within DiffBind, at the cutoff of ≤ 5% FDR and ≥ twofold change. Differential peaks were assigned to the nearest genes using the closestBed command in BEDTools.
Copy number alteration (CNA) analysis. CNA analysis was done according to 25 with minor modifications. Instead of using 1-Mb windows with 500-kb overlap and an absolute z-score cutoff of > = 2, we split each autosome into 100-kb windows with a step size of 10 kb, and accordingly used a more stringent cutoff of absolute z-score > = 3. For each of the normal and T-PLL cases, alignments from H3K4me3, H3K27ac and input libraries were combined and those in the peaks and blacklisted regions (see above) were excluded. BEDTools was used to estimate per window read pairs in each sample. Windows with no (6.37% of the windows) or less than 100 read pairs (0.87% of the windows) were excluded. The coverage for the remaining windows was calculated as: Coverage = (Window read counts scaled to 10 M)/(Effective window size in kb), where effective window size does not include the base pairs overlapping peaks and blacklisted regions. The coverage in each of the six T-PLL was normalized to the average of the median coverage from the three normal, and log2 transformed. Z score was calculated using the above median-normalized coverage in each PLL sample. A copy number gain or loss was inferred if the z score > = 3 or < = − 3, respectively. The differential H3K4me3 and H3K27ac peaks were assigned the CNA status of the nearest 100-kb window with the shortest distance between peak center and window middle point.
TF motif enrichment analysis. For differential H3K4me3 and H3K27ac peaks, we identified the nearest open chromatin regions by intersecting with ATAC-seq peaks from blood 26,27 (see below). Sequences from ATAC-seq peak center ± 50 bp were extracted and enrichment of known TF motifs was identified using the Homer package (v4.10, http:// homer. ucsd. edu/ homer/).

Transcriptional dysregulation in T-PLL.
We performed total RNA sequencing in PBMCs from four T-PLL patients (P1, P3, P5, and P6) and three healthy individuals (N1-N3), and mRNA sequencing for another five T-PLL (P2, P4, and P7-P9) and four normal (N4-N7). Of the nine T-PLL patients, P1 and P3 had already received treatment when samples were collected, while the other seven patients were treatment naive (Table 1). In total RNA sequencing, a total of 12,375 protein-coding genes had RPKM ≥ 0.5 in at least one of the samples. The top 10,000 most variable genes were used in unsupervised clustering, which revealed distinctive differences between T-PLL and normal samples (Fig. 1a). We observed noticeable heterogeneity of gene expression within T-PLL. In mRNA sequencing, a total of 13,728 protein-coding genes had RPKM ≥ 0.5 in at least one of the samples. Unsupervised clustering using the top 10,000 most variable genes revealed a similar pattern of expression changes between T-PLL and normal and heterogeneity across T-PLL patients (Supplemental Fig. 1a).
Using total RNA sequencing, we identified 807 and 865 genes that were up-and down-regulated, respectively, in T-PLL compared to normal (Fig. 1b). In addition, mRNA sequencing revealed 1466 up-and 898 downregulated genes in T-PLL (Supplemental Fig. 1b). Of those, 280 up-and 279 down-regulated genes were shared between the two RNA-seq datasets. mRNA sequencing identified 1.8 × as many up-regulated protein-coding genes in T-PLL compared to total RNA sequencing. This discrepancy is likely due to difference in sequencing depth, platform-specific bias of the two RNA-seq protocols, and the heterogeneity among T-PLL cases. To gain a better understanding of the possible causes, we first estimated the sequencing depth by counting the number of reads mapped to all 57,773 genes and to the subset of 20,327 protein-coding genes. While normal samples had largely comparable number of reads mapped to protein-coding genes in both datasets, T-PLL samples indeed had lower coverage of protein-coding genes in total RNA sequencing (17.87-37.13 M) than in mRNA sequencing (41.52-52.74 M) (Supplemental Fig. 2a). Next, we assessed expression level for the protein-coding genes up-regulated in T-PLL between the two datasets, using gene maximum RPKM values across both normal and T-PLL cases. Up-regulated genes identified in both datasets showed similar levels of expression (Mann-Whitney www.nature.com/scientificreports/ U test p = 0.09, left panel in Supplemental Fig. 2b). Not surprisingly, for the genes identified as up-regulated in T-PLL by mRNA sequencing alone, they had much higher expression levels in mRNA sequencing than in total RNA sequencing (p < 2.2e−16, middle panel in Supplemental Fig. 2b), and vice versa (p < 2.2e−16, right panel in Supplemental Fig. 2b). It is known that genes with low counts are less likely to be significant in differential expression test due to a lack of statistical power 38,39 . For both datasets, Gene Ontology analysis indicated a remarkable enrichment of the down-regulated genes in immune response (Table 2 and Supplemental Table 1). Similarly, pathway enrichment analysis using the Reactome database revealed enrichment in immune and inflammatory responses, including the "chemokine receptors", "immune system", and "immunoregulatory interactions between a lymphoid and a non-lymphoid cell" ( Table 2 and Supplemental Table 1). On the other hand, Gene Ontology analysis revealed that the up-regulated genes were enriched in "WNT-protein binding" and "WNT-activated receptor activity", as well as in many other categories of developmental related biological processes (Supplemental Tables 2 and 3). This result indicates a similar role for WNT signaling network in T-PLL, as documented in B-cell chronic lymphocytic leukemia [40][41][42][43] and T-lineage acute lymphoblastic leukemia 44 .

Global alteration of regulatory regions in T-PLL.
To understand the regulatory landscape and the possible roles of epigenetic mechanisms in T-PLL pathogenesis, we generated ChIP-seq data for H3K4me3 and H3K27ac, which are mainly located in the promoters and active enhancers, respectively. To reveal the extent of global similarities in H3K4me3 and H3K27ac distribution across normal and T-PLL, we estimated pairwise Pearson correlation using normalized read counts from merged peaks present in two or more samples (Supplemental Fig. 3a,b). As expected, the three normal samples were highly similar in the occupancy of both H3K4me3 (R = 0.96-0.97) and H3K27ac (R = 0.90-0.92). On the other hand, the six T-PLL were more divergent (R = 0.82-0.92 for H3K4me3 and R = 0.65-0.86 for H3K27ac), in particular for H3K27ac, suggesting an obvious heterogeneity of epigenomic landscape. Of the six T-PLL, P1 received romidepsin treatment one month before sample collection. Romidepsin is a histone deacetylase (HDAC) class I-selective agent for HDAC1, HDAC2, HDAC3, and HDAC8. Romidepsin treatment is through to increase histone acetylation, thus altering gene expression, in particular activating tumor suppressor genes 25,45 . P1 is a non-responder to romidepsin, and was thus included in the analysis together with the other T-PLL. To support this, we assessed global enrichment of H3K27ac over the www.nature.com/scientificreports/ respective input (Supplemental Fig. 4a), as well as the signal levels within consensus peaks shared either by all six T-PLL or by at least two (Supplemental Fig. 4b). P1 did not show clearly excess enrichment relative to the other T-PLL. Also, we examined the expression of tumor suppressor genes 46 , revealing that P1 had the lowest median expression among the four T-PLL with total RNA sequencing data (Supplemental Fig. 4c).
Specifically, for H3K27ac we focused on distal peaks, defined as those that are ≥ 2.5 kb away from gene TSS; while for H3K4me3, we focused on proximal peaks in the promoter regions, defined as TSS ± 2.0 kb. Using the top 20,000 most variable (out of 52,248) distal H3K27ac peaks that are present in two or more samples, unsupervised clustering identified a preferential loss of active enhancers in T-PLL (Fig. 2a). Similar results were obtained with the top 10,000 and 30,000 distal H3K27ac peaks (Supplemental Fig. 5a,b). There are 16,629 proximal H3K4me3 peaks within TSS ± 2.0 kb. A similar pattern was identified for the 10,000 most variable H3K4me3 peaks (Fig. 2b). Therefore, there is a genome-wide change of regulatory landscape in T-PLL. As expected, we observed obvious heterogeneity in the occupancy of both H3K4me3 and H3K27ac, showing two sub-groups with P2, P4 and P5 belonging to one and the other three forming the second one (Fig. 2a,b and Supplemental Fig. 5a,b).
We next used the DiffBind package to identify the changes of H3K4me3 and H3K27ac between the four T-PLL (P1, P3 and P5-P6) and three normal that also had total RNA sequencing data. The analysis identified differential occupancy for 15.1% (9362) of the 62,016 H3K27ac peaks and 10.5% (3935) of the 37,306 H3K4me3 peaks ( Table 3). As revealed by unsupervised clustering (Fig. 2a,b), there is a strong tendency of losing both marks in T-PLL compared to normal. To understand whether differential peaks are associated with the change of gene expression, we assigned them to the nearest genes and checked the expression changes of those genes. Indeed, for both marks, peaks that showed increased occupancy in T-PLL were preferentially associated with up-regulated in T-PLL, compared to non-differential peaks, and vice versa (enrichment: 3.8-19.2, chi-square test p < 2.2e−16) (Table 3). To correlate the changes of gene expression with those of proximal occupancy of H3K4me3 and H3K27ac, for the 865 down-and 807 up-regulated genes identified by total RNA sequencing, we also plotted the average signal profiles of H3K4me3 and H3K27ac over the TSS ± 2 kb (Fig. 2c,d). As expected, overall, the down-regulated genes were associated with decreased H3K4me3 and H3K27ac, and vice versa; the correlated changes were more obvious for the up-regulated genes (Fig. 2d).
To understand which TFs might be involved in the perturbation of T-PLL epigenome, we first used diffbind to identify differential H3K4me3 and H3K27ac peaks on chr1-22 between the six T-PLL and three normal, followed by k-means clustering of the differential peaks. The analysis identified 2445 (1728 + 717) of the 2456 H3K27ac peaks with increased signal in T-PLL and all 5587 peaks with decreased signal as discriminatory peaks (Fig. 3a). To a less extent, 80% (647 + 124 + 129 = 900) of the 1128 H3K4me3 peaks with increased signal in T-PLL and 90% (2658) of the 2971 peaks with decreased signal were retained (Fig. 3b). The other differential peaks showed less difference between T-PLL and normal. We then used the Homer package to identify known TF motifs that are enriched in these discriminatory peaks (Fig. 3c,d). The differential H3K27ac peaks (Fig. 3c), as well as H3K4me3 peaks gained in T-PLL (upper panel in Fig. 3d) are enriched for binding motifs of E-twenty-six (ETS) TFs, such as ERG, ETS1, and ELF4. On the other hand, differential H3K4me3 peaks lost in T-PLL (lower panel in Fig. 3d) are particularly enriched for Interferon regulatory factor (IRF) TF motifs. Finally, for both marks, the differential Table 2. Enriched pathways and GO terms for genes down-regulated in T-PLL. Gene expression was quantified by total RNA sequencing for three normal (N1-N3) and four T-PLL (P1, P3, P5, and P6). www.nature.com/scientificreports/ To reveal the extent by which CNA may contribute to epigenomic change, we identified CNAs in individual T-PLL cases based on sequence coverage in 100-kb windows along chr1-22. At the cutoff of absolute z-score > = 3, there are on average 0.34% and 0.37% of the windows that showed CN gain and loss, respectively. However, peaks are > sevenfold over-represented in the windows with CN gain compared to those with CN loss. Specifically, for the consensus peaks shared by at least two T-PLL, 0.78-5.20% of the H3K27ac peaks and similar proportions of the H3K4me3 peaks (0.66-4.82%) fall in windows with CN gain, versus < = 0.6% in windows with CN loss (Supplemental Table 4). This bias is unlikely due to peak calling, as a similar trend was also observed for additional two consensus peak sets of each mark, i.e., those present in > = 2 T-PLL and > = 2 normal, as well as those present in > = 2 normal but absent from all T-PLL. In these additional peak sets, < 0.5% of peaks overlap windows with CN loss, changing little (< 1%) at a less stringent cutoff of z-score < = − 2. The results indicate that regions with CN loss are less likely to be enriched with H3K4me3 and H3K27ac modifications. Finally, we investigated the overlap of differential peaks with windows carrying CNAs (Supplemental Table 4). Intersecting both marks with windows carrying CN gain revealed overlap with about 3-4% the peaks with increased signal in T-PLL, versus with 0.6% (H3K27ac) and 1.94% (H3K4me3) of the peaks with decreased signal in T-PLL. Thus, even for peaks falling within windows with CN gain, some showed reduced histone modification levels in T-PLL instead. On the other hand, similar to the whole peak set, differential peaks are much less associated with windows with CN loss. It is possible that, with a large (100-kb) window size and a modest coverage, small CNAs overlapping differential peaks may not have been identified, and that this method may be less effective in the detection of CN loss compared to that of CN gain. Nevertheless, CNAs were found to contribute to the changes of H3K4me3 and H3K27ac in a small portion of peaks, including those associated with key genes (such as MYC and AGO2) described below.
Changes in chromatin states of key dysregulated genes. As revealed by a recent study 5 and our RNA-seq data described above, T-PLL is characterized by the dysregulation of several oncogenes, such as TCL1A, MYC and EZH2, as well as genes implicated in DNA damage repair (e.g., ATM) and T-cell receptor regulation (such as CTLA4 and SLAM family members). Below we focused on these key genes. Specifically, we aimed to illustrate the epigenetic state in the regulatory regions of these genes and to understand how epigenetic alterations have modulated the transcriptional outcome.
The hallmark of genetic alterations in T-PLL includes t(14;14)(q11;q32) and inv (14)(q11q32) that involve the proto-oncogene TCL1A. TCL1A activates the serine/threonine kinase AKT that is implicated in cell survival regulation, as well as the NF-κB signaling pathway that plays key roles in the initiation and promotion of cancer 6,50 . Increased TCL1A expression is associated with less favorable clinical outcome in chronic lymphocytic leukemia 6,51 . TCL1A was significantly up-regulated in T-PLL compared to normal (log 2 FC = 5.60 and p = 3.30e−08 in total RNA-seq; log 2 FC = 7.87 and p = 1.08e−18 in mRNA-seq) (Fig. 4a, right panel). All six T-PLL cases demonstrated a gain of H3K4me3 and H3K27ac over its promoter and gene-body (Fig. 4a), which, we believe, contributes to TCL1A activation in addition to the known rearrangement of TCL1A locus. Based on the signal level in the corresponding input, at least two of the T-PLL samples (P1 and P5) showed copy number gains in TCL1A, which may have also resulted in increase of TCL1A expression, as well as of H3K4me3 and H3K27ac occupancy. Table 3. Association of differential H3K27ac and H3K4me3 occupancy with differential gene expression. Total RNA-seq differential expression analysis and ChIP-seq DiffBind were performed on four T-PLL (P1, P3, P5 and P6) and three normal (N1-N3). Peaks were split into those that showed differential occupancy between T-PLL and normal and others that showed no differential occupancy based on DiffBind analysis (FDR ≤ 0.05, fold-change ≥ 2). Peaks were assigned to the nearest genes. For differential peaks showing increased occupancy in normal ("Up in normal"), the number in parentheses indicates the number of up-regulated genes in normal. For differential peaks showing increased occupancy in T-PLL ("Up in T-PLL"), the number in parentheses indicates the number of up-regulated genes in T-PLL. Enrichment level is estimated as: (Number of differential peaks linked to differentially-expressed genes/total differential peaks) / (Number of non-differential peaks linked to differentially-expressed genes/total non-differential peaks).

Mark Region
Differential peaks www.nature.com/scientificreports/ www.nature.com/scientificreports/ www.nature.com/scientificreports/ To examine whether this ~ 6-kb H3K27ac peak exists in other tissue and cell types, we first checked the 98 human reference epigenomes that have H3K27ac data 28 . H3K27ac enrichment was detected only in a primary B-cell sample and in GM12878, which is a B cell-derived lymphoblastoid cell line. Based on other publicly available H3K27ac ChIP-seq data we collected, this broad peak was identified in four of the five normal B-cell samples, in patient biopsies from all seven mantle cell lymphoma, three small lymphocytic lymphoma and four of the seven high-grade B-cell lymphoma subjects (Supplemental Fig. 6), but not in any of the 10 adult T-cell leukemia/ lymphoma (ATL) patients 35 . In addition, based on the catalog of super-enhancers from 86 cell and tissue types 36 , this enhancer is part of a B-cell specific super-enhancer. These results suggest a possibility that requisition of this cell-type restricted enhancer activates the oncogenic expression of TCL1A in T-PLL.
Oncogenic activation of MYC in T cell leukemia 52 and B cell lymphoma 33 involves interaction between its promoter with subtype specific enhancers. T-PLL is also characterized by the activation of MYC 5 . We detected upregulation of MYC expression (log 2 FC = 2.83 and p = 5.13e−03 in total RNA-seq; log 2 FC = 3.37 and p = 9.49e−06 in mRNA-seq). To support a role for epigenetic mechanism in MYC activation, ChIP-seq data showed an overall increase of H3K4me3 and H3K27ac occupancy over an ~ 6 kb region in all T-PLL cases that covers MYC promoter and gene body (Fig. 4b). The signal in the inputs suggested copy number gains in P1, P5 and P6. To identify the enhancer that potentially interacts with the MYC promoter, we used the promoter CHi-C data from 17 blood cell types 37 . The analysis identified an ~ 50-kb enhancer located over 500 kb upstream that interacted with the MYC promoter in total B, naïve B, fetal thymus and five of the six types of T cells (Fig. 4b). This broad enhancer showed marked increase of H3K27ac in three of the T-PLL cases (P1-P3). Thus, in these three cases, the up-regulation of MYC expression appears to be driven by the increased H3K27ac in the enhancer together with increased H3K4me3 and H3K27ac in the promoter. For the other three cases (P4-P6), the increased MYC expression is associated with a gain of H3K27ac in the promoter in all three, and a gain of H3K4me3 in P6.
EZH2 is an oncogene encoding a histone methyltransferase that deposits H3K27me3, a histone mark associated with transcriptional repression 53 . EZH2 is highly expressed in several B-cell lymphomas 53 . We found that EZH2 was up-regulated in four of the T-PLL samples (P2-P4 and P6) (Supplemental Fig. 7). Three of the four (P3, P4, and P6) also showed increased occupancy of both H3K4me3 and H3K27ac in the promoter (Supplemental Fig. 7), with the exception of P2 that had a copy-number loss indicated by the lower signal in the input. Overall, the elevated EZH2 expression in T-PLL is associated with increased H3K4me3 and H3K27ac in its promoter.
In the argonaute (AGO) family, only AGO2 has the catalytic activity and plays an essential role in small RNAguided gene silencing 54 . Often up-regulated in cancer, AGO2 is involved in chromatin remodeling and alternative splicing 5 . AGO2 was up-regulated in P1 and P2 (Supplemental Fig. 8). All the six T-PLL samples showed increased H3K4me3 and H3K27ac in its promoter (Supplemental Fig. 8). The signal level in inputs suggested that four (P1, P3, P5 and P6) had copy-number gains, which likely have contributed to the increased H3K4me3 and H3K27ac in these cases. For P2 that showed no copy number alteration, its up-regulation may be directly due to the increased H3K4me3 and H3K27ac.
On the other hand, we also checked the chromatin states for several T-cell receptor regulators and ATM. As a major negative regulator of T-cell response, CTLA4 was markedly down-regulated in T-PLL (log 2 FC = − 6.96 and p = 8.46e−14 in total RNA-seq; log 2 FC = − 6.27 and p = 1.08e−16 in mRNA-seq), which was associated with an obvious reduction of H3K4me3 and H3K27ac (Fig. 4c). Input signal indicated no CNV in any of the T-PLL cases. The SLAM family members also act as the negative regulators of T-cell response. Of the nine members, seven reside in a 378-kb region on chromosome 1 (Supplemental Fig. 9). Of the seven, CD84 (also known as SLAMF5, log 2 FC = − 3.69 and p = 2.22e−05 in total RNA-seq; log 2 FC = − 4.47 and p = 2.85e−10 in mRNA-seq), SLAMF7 (log 2 FC = − 3.12 and p = 4.44e−04 in total RNA-seq; log 2 FC = − 3.78 and p = 7.42e−08 in mRNA-seq) and CD244 (SLAMF4, log 2 FC = − 2.17 and p = 2.26e−02 in total RNA-seq; log 2 FC = − 2.34 and p = 1.16e−03 in mRNA-seq) were down-regulated in all T-PLL cases. Accordingly, these three genes showed lower levels of H3K4me3 in their promoters compared to the normal. In addition, SLAMF6 was down-regulated in P2, P4 and P5, showing lower H3K4me3 in promoter in P4 and P5; SLAMF1 was down-regulated in P1, P3 and P4, with lower H3K4me3 in promoter in the first two. Clearly, for genes involved in the regulation of T-cell response, alterations of chromatin state in the promoters are linked to their dysregulation.
In B-cell chronic lymphocytic leukemia, ATM was found to physically interact with TCL1A protein, and this association activates the NF-κB pathway in the regulation of apoptosis 6,55 . ATM aberration is viewed as a key initiating lesion in T-PLL 5 . ATM showed lower H3K4me3 occupancy in the promoter in all T-PLL cases, with down-regulation of expression in four of them (except in P1 and P5) (Supplemental Fig. 10). Based on the input signal within ATM, there was a copy-number loss in five T-PLL cases (P2-P6), which may have led to the observed changes in epigenetic state and expression. The result suggested a link between the reduction of H3K4me3 in promoter and ATM down-regulation. Collectively, our analysis demonstrated a central role for epigenetic mechanisms in the dysregulation of key oncogenes, and genes with roles in DNA damage repair and T-cell receptor regulation in T-PLL.

Discussion
T-PLL was best characterized by oncogenic expression of TCL1A and dysfunction of ATM 5 . Our analysis validated the previous genomic findings of dysregulated genes and further suggested key roles for epigenetic mechanisms in transcriptional dysregulation. We identified genome-wide alterations of chromatin states at both promoters and active enhancers in T-PLL, which likely involves immune related TFs that are highly enriched in differential peaks. The results supported the epigenetic regulation of oncogene activation and repression of genes involved in DNA damage and T-cell responses. Importantly, we provided evidence for the gained enhancers in T-PLL that contribute to the overexpression of TCL1A and MYC. We believe the same mechanism underlies their activation in several other hematological malignancies. For example, a super-enhancer within TCL1A was www.nature.com/scientificreports/ identified in high-grade B-cell lymphoma, mantle cell lymphoma and small lymphocytic lymphoma 33,34 . Accordingly, TCL1A showed overexpression in these B-cell lymphomas 56,57 . In B-cell chronic lymphocytic leukemia, TCL1A up-regulation was also associated with the gained super-enhancer 58 and with DNA hypomethylation in its promoter 55 . Additionally, in many other cancers including T cell leukemia, chronic myeloid leukemia and multiple myeloma, MYC expression was linked to the super-enhancers acquired in tumor cells 36 . Thus, this study has enhanced our understanding of epigenetic regulation of gene expression in T-PLL. Our study is limited by a small sample size because T-PLL is a rare mature T cell leukemia. Nevertheless, we revealed an obvious heterogeneity in gene regulatory landscape within T-PLL, with P2, P4 and P5 in one group and the other three in the second group. We also identified marked differences in gene expression across T-PLL cases. Further study using a larger collection of samples is needed to identify the epigenetic and transcriptional programs associated with disease progression.
Our analysis supported roles for epigenetic mechanisms in gene dysregulation in T-PLL. Genetic alterations were found to contribute to transcriptional abnormality in T-PLL 5 . Based on input signals, two to four of the T-PLL cases had copy-number gains in AGO2, MYC and TCL1A that were accompanied by increased H3K27ac and H3K4me3, while copy-number loss was observed in ATM in five T-PLL cases (P2-P6) that was accompanied by decreased H3K4me3. Complementarily, CNA analysis with 100-kb windows identified CN gains covering AGO2 and MYC, and also TCL1A at a less stringent cutoff of z = 2. For these four genes, the CNVs inferred from input signal spanned both gene-body and promoter. In fact, copy-number gain in AGO2 and MYC and copynumber loss in ATM are the major genomic aberrations in T-PLL 5 . Genomic rearrangements in regulatory regions represent an important mechanism for oncogenic activation. For example, in acute myeloid leukemia 59 and acute lymphoblastic leukemia 60 , the distal super-enhancers, 1.7 and 1.3 Mb away from the MYC promoter, respectively, were located in regions of focal amplification. Integrative analysis of genetic variation, epigenetic profiles with gene expression will help decipher the roles of interplay between genetic and epigenetic factors in shaping the transcriptional program in T-PLL.
In summary, the analyses have identified the transcriptional changes for genes enriched in immune and inflammatory responses, as well as for those in WNT signaling network. Using H3K4me3 and H3K27ac ChIPseq data, we revealed genome-wide alterations of chromatin states at both promoters and active enhancers in T-PLL. Integrative analysis with RNA-seq data suggests the changes of regulatory regions as a major mechanism underlying transcriptional reprogramming. Importantly, epigenetic changes in regulatory regions are evident at key oncogenes and genes involved in DNA damage response and T-cell receptor regulation, which are tightly linked to their transcriptional dysregulation, the core lesions of T-PLL. Collectively, this study enhances our understanding of the roles of epigenetic mechanisms in T-PLL pathogenesis.