
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)12530-3
10.1016/j.heliyon.2024.e36499
e36499
Research Article
PML is a constitutive component of chromatin domains enriched in repetitive elements and duplicated gene clusters in cancer cells
Fracassi Cristina ab
Simoni Matilde a
Uggè Martina a
Morelli Marco J. b
Bernardi Rosa bernardi.rosa@hsr.it
a⁎
a Division of Experimental Oncology, IRCCS San Raffaele Scientific Institute, Milano, Italy
b Center for Omics Sciences, IRCCS San Raffaele Scientific Institute, Milano, Italy
⁎ Corresponding author. bernardi.rosa@hsr.it
17 8 2024
15 9 2024
17 8 2024
10 17 e3649928 3 2024
10 8 2024
16 8 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).
Heterochromatin is a pivotal element in the functional organization of genomes. In our study, we delve into the heterochromatin pattern of association by the PML (promyelocytic leukemia) protein. By using PML chromatin immunoprecipitation and sequencing data and comparing computational methodologies to depict PML chromatin association, we describe PML-associated domains or PADs as large heterochromatic regions that exhibit similar genomic features across cancer cell lines. We show that PADs are specifically enriched in non-coding genes, duplicated gene clusters, and repetitive DNA elements. Moreover, we find enriched binding motifs of KZFPs, which are involved in orchestrating epigenetic repression at repetitive DNA elements. Hence, our findings suggest that PML conservatively associates to heterochromatic domains enriched in repetitive DNA elements and duplicated gene clusters in cancer. These findings contribute to a broader understanding of the complex regulatory framework of genome organization by heterochromatin in cancer.

Highlights

• Computational analysis of PML ChIP-seq data relies on the identification of broad DNA association.

• Repetitive DNA elements are hallmarks of PML association to DNA in human cancer cells.

• PML emerges as a key component of similar chromatin domains enriched in duplicated gene clusters and repetitive DNA elements.
==== Body
pmc1 Introduction

The PML (promyelocytic leukemia) nuclear protein provides a molecular scaffold for the nucleation of insoluble biomolecular condensates named PML nuclear bodies (PML-NBs) [1]. PML-NBs are generally described as platforms for the regulation of proteins therein associated upon cellular stress. Several transcription factors and chromatin regulators localize to or transit through the PML-NBs, denoting an important role of PML in the regulation of transcription and genome functions [2]. Accordingly, although the PML protein does not contain DNA-binding motifs, it has been variously associated to transcribed DNA, acetylated blocks of chromatin and specific genomic loci via antibody-based imaging studies that focused on PML-enriched nuclear bodies [2]. This prompted researchers to characterize chromatin regions in contact with the PML-NBs at the whole genome level. Chromatin immunoprecipitation (ChIP) is one of the most widely used techniques for the identification of DNA sequences associated to a given protein. However, ChIP optimally works with soluble proteins while the PML-NBs are highly insoluble aggregates and are underrepresented in ChIP fractions. Moreover, because PML also exists in a nucleoplasmic or non NB-bound conformation [1], since ChIP an antibody-based approach, it cannot discriminate between the nucleoplasmic pool or the PML-NB-bound PML fraction [2]. To overcome this limitation several techniques have been designed to enrich DNA associated to the PML-NBs via proximity labeling [[3], [4], [5]]. These approaches (ALAP-seq and Immuno-TRAP) leverage on the aggregation of PML moieties within the PML-NBs to label chromatin regions that are proximal to these structures and led to the identification of narrow peaks of promoter-enriched genomic regions at specific loci, like the p53 and MHC genes [2]. These studies suggested that the PML-NBs exert functions of transcriptional regulation in nearby genomic regions.

In this scenario, the genomic characterization of non NB-bound PML remains poorly described and only few studies have reported an association of PML to DNA outside the PML-NBs via ChIP-seq and DNA FISH approaches [6,7]. Of note, it was reported that common peak calling algorithms designed for transcription factors could not identify regions of PML chromatin association from ChIP-seq experiments [5], suggesting that the fraction of PML that is immunoprecipitated with this biochemical approach associates to DNA differently. Accordingly, by performing ChIP-seq experiments in mouse and human cells we and others have positioned PML to large heterochromatic domains named PML-associated domains or PADs, in a DNA association profile that is broad and of low intensity [6,7]. These data imply that PML connects to chromatin in a manner akin to chromatin structural proteins like lamins. DNA FISH coupled to visualization of PML-NBs by immunofluorescence showed that PADs are at large distant from PML-NBs [6,7]. However, a rigorous computational analysis of PML ChIP-seq experiments within the same cellular context is still lacking and no standardized way of analyzing PML chromatin association exists.

In this work, we dissected the profile of PML DNA association in ChIP-seq experiments from different cancer cell lines. We benchmarked different peak calling algorithms and describe a methodological approach to identify chromatin regions of reliable PML association on ChIP-seq data. Through a comparative analysis of PML genome-wide binding profiles in cancer cell lines, we show that PML associates to broad heterochromatic domains that are enriched in repetitive elements and functionally related duplicated gene clusters, suggesting that PML may participate to the regulation of chromatin states at specific genomic regions.

2 Materials and methods

2.1 Cell culture, lentiviral vectors, lentiviral production and transduction

MDA-MB-231 and RCC4 cells were purchased from ATCC and maintained in DMEM supplemented with 10 % fetal bovine serum (FBS) (Euroclone) and 1 % Penicillin/Streptomycin antibiotics (Lonza) at 37 °C in a humidified atmosphere containing 5 % CO2. Third generation lentivirus (LV) stocks were prepared, concentrated and titrated as previously described [26]. Briefly, self-inactivating (SIN) LV vectors were produced by transient transfection of HEK293T cells with the packaging plasmid pMDLg/pRRE, Rev-expressing pCMV-Rev, the VSV-G envelop-encoding pMD2.VSV-G plasmids, and shRNA-carrying vectors. Optimal puromycin concentration was pre-determined by performing dose-response curves and used at a final concentration of 2.5 μg/ml for MDA-MB-231 and 1 μg/ml for RCC4. In RCC4 cells. In RCC4 cells, inducible silencing of PML was induced with doxycycline monohydrate (Merck, D1822-500 MG) and at a final concentration of 100 ng/ml for 96 h. Culture media was changed every 48 h.

2.2 Chromatin-immunoprecipiation sequencing (ChIP-seq)

Cells were seeded 24 h before chromatin isolation at 60–70 % confluence in 15 cm plates. Cells were trypsinized and double crosslinking was performed in suspension. Cells were first resuspended in PBS containing 2 mM Di(N-succinimidyl) glutarate (DSG, Sigma-Aldrich 80424) for 45 min at RT on gentle rotation. Cells were then centrifuged at 1200 rpm at RT for 5 min. The resulting pellet was resuspended in PBS containing 1 % Formaldehyde (Sigma-Aldrich 252549) for 10 min at RT on gentle rotation. After formaldehyde quenching with Glycine (final concentration 125 mM) for 5 min, cells were centrifuged at 1350×g at 4 °C for 5 min, and the supernatant was discarded. Chromatin extraction was performed as previously described [27] and sonication was performed using the Bioruptor (Diagenode Bioruptor 300) at high intensity, 30 s ON and 40 s OFF for 12 cycles for MDA-MB-231 cells and 30 s ON and 50 s OFF for 17 cycles for RCC4 cells, to obtain chromatin enriched in fragments of 200–1500 bp. For chromatin quantification, an aliquot of chromatin was de-crosslinked and purified using QIAquick PCR purification kit (Qiagen). Chromatin was then quantified with Nanodrop spectrophotometer and efficiency of sonication was measured by agarose gel electrophoresis. For sequencing, DNA quality was evaluated with a High Sensitivity D5000 ScreenTape (Agilent Technologies). After sedimentation, 50 μg (for sequencing) of chromatin were incubated with 20 μg of anti-PML (Santa Cruz Biotechnology 71910) and 5 % of chromatin was collected as input sample. Chromatin was pre-coupled to target-antibody overnight at 4 °C on rotation and the day after beads were added (ChIP-IT Protein G magnetic beads 53033 or Magna ChIP™ Protein A Magnetic Beads 16–661) for 4 h at 4 °C on rotation. ChIP samples were washed three times in ice-cold RIPA buffer +1 % SDS (Tris–HCl pH 8 10 mM, NaCl 140 mM, Triton X-100 1 %, Na-deoxycholate 0.1 %, EDT A 1 mM, EGT A 0.5 mM) and elution was performed in 250 μl of Elution buffer (NaCl 50 mM, Tris–HCl pH 7.5 20 mM, EDTA 5 mM, SDS 1 %, RNAse-A 0.5 μg/ml, Proteinase K 2 μg/ml) for 6 h at 37 °C on rotation. Finally, beads were removed, and samples were de-crosslinked overnight at 65 °C. DNA was purified with the QIAquick PCR purification kit (Qiagen) and used to construct libraries following with the ChIPSeq Illumina protocol (Illumina). After being barcoded, pooled and sequenced libraries where on sequenced an Illumina Nova-Seq 6000 system. ChIP-seq experiments were performed generating 40 million reads, 100 nucleotide long, in paired end. Of note, while 40 million sequencing reads were selected to be sufficient to examine PML enrichment in ChIP-seq experiments, we are aware that increasing the number of reads may highlight additional regions of PML enrichment. After sequencing, reads were trimmed using BB-Duk from BBTools suite version 37.36 with suggested settings (ktrim = r k = 23 mink = 11 hdist = 1) and then mapped using BWA-MEM version 0.7.17 [28] on the human genome assembly GRCh38. Uniquely mapped reads were selected with markDuplicates from Picard Tools version 3.1.1.

2.3 Peak calling

To map PADs, ChIP-seq normalization bias were avoided by down-samplings for each chromosome each pair of mapped ChIP and input read files. Mapped reads were used to call domain using ten runs of Enriched Domain Detector (EDD) [8] with auto-estimation of GapPenalty and BinSize, and mean GapPenalty and BinSize values from these runs were used for a last run. Final domains were the union of domains of all replicates. To identify PML associated peaks, peak calling was performed with MACS2 version 2.2.6 [10] (with BAMPE setting and qvalue cutoff = 5.00e-02), SICER version 1.0.3 [9] (with window size = 200, gaps size = 600, fragment size = 300 and fdr cutoff = 0.1), Peakranger version 1.18 [11] (with default settings). Further filtering was done on peaks mapping in regions present in the ENCODE hg38 blacklist [29].

2.4 Repeat Enrichment Estimator

PML ChIP-seq and matched input samples were aligned to a repeat database with Repeat Enrichment Estimator (RepEnrich) [30]. The annotation files used were repeatmasker files with simple and low-complexity repeats removed (with satellite repeats and transposons included). To identify uniquely and multimapping reads, data were mapped with map Bowtie version 1.3.1 [31] and counts of repetitive elements were calculated by RepEnrich [30] using default settings. The final output estimating counts of repeats (fraction of counts) was then used to build a table of counts for all conditions (PML ChIP and Input samples) and differential expression analysis was performed with DeSeq2 package [32] to identify differential enriched repetitive elements in PML ChIP-seq compared to input samples. Specifically, DeSeq2 analysis was used to obtain a log2 fold change values for PML ChIP-seq with respect to input and an associated p-values for each repetitive element. ggplot2 [33] was used for plots.

2.5 RNA-sequencing (RNA-seq)

RNA sequencing experiments were previously described [[7], [12]], generating 30 million single end reads, 100 nucleotide long. Sequencing adapters were removed using trimmomatic v0.39 [34], and fastq files were then aligned to the human genome assembly GRCh38 (hg38) using the STAR aligner v2.5.3a [35]. Annotation of genomic features was performed using the feature-Counts tool v1.6.4 [36] using the GENCODE v31 Gene transfer format (GTF). Differential gene expression was evaluated in R/BioConductor using the DeSeq2 package [32] and using a false discovery rate (FDR) of 0.05 for significance. Comparisons were performed between cells expressing PML-specific shRNA (shPML) and control shRNA (shCTRL).

2.6 Operations on genomic intervals

Processing of peaks and domains was performed using BEDTools v2.30 [37] and BEDOPS v2.4.41 [38]. The number of overlapping peaks/domains between different conditions was computed with the intervene venn function from the Intervene v0.5.8 package [39]. Common and unique peaks/domains were identified with bedops intersect or bedtools substract. Mean Log2(ChIP/Input) or normalized counts were calculated with bamCompare (-bs 1000 --scaleFactorsMethod SES) or bamCoverage (-bs 50 --normalizeUsing RPGC) and quantified using multiBigwigSummary from Deeptools version 3.5.1 [40]. Otherwise, bigwig tracks were generated from Log2(Chip/Input) ratios in 1-kb bins using EDD [8]. All genomic intervals and profiles were visualized using Integrative Genomics Viewer [41]. To assess ChIP data quality and reproducibility, Pearson correlations were determined between peaks from each replicates with the intervene pairwise function from Intervene v0.5.8 package [39]. Genes were ascribed to an EDD called domain if they overlapped with a domain by at least one base-pair.

2.7 HOMER motif analysis

For transcription factors binding analysis of genes regulated by PML in PADs, motif calling was performed by considering a 10 kb region upstream the TSS of the gene of interest. Common PADs were fragmented in 1 kb sized regions with the bedops chop option from BEDOPS v2.4.41 [38]. HOMER version 4.11 [42] was run with the following manner: findMotifsGenome.pl input.txt hg38 -size given. We considered as enriched motifs results from the HOMER de novo motif discovery step (homerResults).

2.8 Functional enrichment analysis

For RNA-seq data, differential gene expression sets were ranked based on their fold change values and subsequently analyzed with fgsea version 1.28.0 [43]. In addition, gene set enrichment analysis was performed with the ShinyGO 0.80 webtool [44] and KEGG was used as reference databases to identify significant pathways with FDR cutoff <0.05.

2.9 Karyotype analysis

Distribution of genes falling within common PADs was visualized by using karyoploteR [45]. Evaluation of significantly enriched gene clustered was performed by ClusterLocator [46], which perform a two-sided Kolmogorov-Smirnov test on each chromosomal segments to test if the analyzed genes are uniformly distributed. In addition, comparison is done between the input data and with the clustering found in 1000 lists of random genes sets.

2.10 Statistical analysis

Data were processed using GraphPad Prism version 9.0.2 (GraphPad Software, San Diego, California, USA, www.graphpad.com), and the R statistical environment.

3 Results

3.1 PML associates to broad chromatin domains

ChIP-seq data analysis is an inherently complex process due to the underlying variety of DNA association specificities by different proteins. This has led to the development of various peak calling algorithms with distinct statistical frameworks and computational methodologies. As a poignant example, the broad genomic association of chromatin structural proteins such as lamins differs from that of transcription factors (TFs), which bind DNA at restricted gene regulatory regions. Hence, different peak calling algorithms have been devised to capture such diversity. Peak calling algorithms that are used to identify sharp regions of high DNA binding intensity, such as that of TFs, include MACS (Model-based Analysis of ChIP-seq) [10] and PeakRanger [11], while algorithms like SICER (Spatial Clustering for Identification of ChIP-Enriched Regions) [9] and EDD (Enriched Domain Detector) [8] identify broader and low-level enrichment that is representative of DNA association by structural proteins. In addition, EDD was designed to focus on width of enriched genomic regions rather than enrichment strength and provides increased robustness against local variations [8]. By using EDD, we and others have identified broad domains of PML chromatin association from ChIP-seq experiments in mouse and human cells [6,7]. However, a rigorous comparative analysis of PML chromatin association by using different peak calling algorithms is still lacking.

With this in mind, we compared different peak calling algorithms on published PML ChIP-seq data and matched input sequences obtained from MDA-MB-231 triple-negative breast cancer (TNBC) cells and deposited at NCBI GSE226060 [7]. Specifically, we compared the output of two broad (EDD and SICER) and two narrow (MACS2 and PeakRanger) peak calling algorithms. For simplicity, in the description of our analysis we refer to genomic regions identified by all algorithms as ‘peaks’, even though some of these regions are large domains rather than narrow and sharp peaks.

Total genome coverage under peaks detected by these algorithms varied from an average of ∼68–900 Mb (SICER and EDD respectively) to ∼0,2–5 Mb (MACS2 and PeakRanger respectively; Table S1). Importantly, while peaks identified with EDD were conserved between replicates, peaks overlaps were sizably lower with the other peak calling algorithms, especially MACS2 and PeakRanger (Fig. 1A and Fig. S1A). Consistently, Pearson correlation analysis revealed that only EDD peaks have a correlation coefficient >0.6 (Fig. S1B), also upon cross correlation with peaks called by the other algorithms (Fig. S1C), suggesting that EDD is the most reliable method for identifying PML-associated domains. Accordingly, inspection of PML ChIP-seq profiles revealed a conserved broad and low intensity enrichment of PML over input that can be visualized as ratios of ChIP/input reads when measured in 1 kb bins throughout the genome (Fig. 1A and Fig. S1D). In contrast, the other peak calling algorithms, which are designed to detect regions of high enrichment in smaller chromatin regions, did not identify conserved peaks among replicates (Fig. S1D).Fig. 1 Comparison of peak calling algorithms for the analysis of PML ChIP-seq in MDA-MB-231 cells. (A) Genome browser view of peaks identified with the indicated peak calling algorithms in all PML ChIP-seq replicates. (B) Enrichment of PML in EDD, SICER, MACS2 and PeakRanger identified peaks; bar, median; whiskers, min-max; ***P < 10−3. One-way analysis of variance, Dunnett's multiple comparison test. (C) FRIP within EDD, SICER, MACS2 and PeakRanger identified peaks. (D) Venn diagram of overlapping peaks between EDD, SICER, MACS2 and PeakRanger identified peaks. (E) Enrichment of H3K9me3 and H3K27me3 (left graph) and H3K4me3 and H3K27ac (right graph) in PML associated peaks; bar, median; whiskers, min-max; ***P < 10−3. One-way analysis of variance, Dunnett's multiple comparison test.

Fig. 1

To compare regions of PML association identified by different peak callers, we proceeded with the analysis of peaks that were conserved within triplicates (intersecting peaks), even though for MACS2 and PeakRanger these represented a small number of peaks (Fig. S1A). We measured PML enrichment within peaks and found that peaks identified by EDD show the highest enrichment of PML over input, with all enrichment values above 0 (Fig. 1B). In contrast, peaks identified by SICER, MACS2 and PeakRanger display a wide range of enrichment levels, with an average enrichment over input at or below 0 (Fig. 1B). Accordingly, inspection of fraction of reads in peaks (FRIP) indicated that amongst all identified peaks, only EDD had libraries with low background signal and overall concordance between replicates (Fig. 1C). Also, the few common peaks identified by MACS2 within triplicates map to regions of non-specific high read signal (similar to input) that are at large not contained within canonical chromosomes of the hg38 genome, but map to random chromosomes (Table S2) and contain PML signals that are on average equal or lower than the input (Figs. S2A and B). Similarly, peaks identified by PeakRanger, albeit mapping to hg38 canonical chromosomes (Fig. S2A and Table S2), also mark regions where PML signals are overall lower than the input (Fig. S2C). These data suggest that most peaks identified by SICER, MACS2 and PeakRanger are false positive. Consistently, the majority of peaks identified by MACS2 and PeakRanger did not overlap with EDD domains, and only 27 % of peaks identified by SICER overlapped with EDD domains (Fig. 1A and D).

Finally, we used histone modifications to describe the qualitative nature of PML-associated chromatin and confirmed that EDD identified heterochromatic regions enriched in H3K9me3 and H3K27me3, as previously observed [6,7], while PML peaks identified by SICER, MACS2 and PeakRanger were enriched in euchromatic histone modifications H3K4me3 and H3K27ac (Fig. 1E).

Thus, by comparing different peak calling algorithms we show that narrow peak calling is not suitable to identify PML chromatin association in ChIP-seq experiments and we confirm PML association to broad heterochromatic domains [6,7]. In concluding, it is worth mentioning that although PML ChIP-seq protocols do not differentiate between the PML-NBs and the nucleoplasmic, free PML pool, most likely they enrich for the PML free pool because clearing of the cell lysates will presumably deplete the insoluble pool of PML aggregated in PML-NBs.

3.2 PADs identify chromatin domains with similar genomic features

Association of PML to large heterochromatic domains named PADs was first described in mouse embryonic fibroblasts (MEFs) [6] and more recently confirmed in TNBC cells, where PML displays oncogenic functions [7]. To understand if the genomic features of PADs are comparable in other cancer cells where PML is overexpressed and plays oncogenic functions, we mapped PML chromatin association in RCC4 cells, a cell line representative of clear cell renal cell carcinoma (ccRCC) where PML drives cell proliferation and tumor progression [12].

EDD identified megabase-sized domains also in RCC4 cells (Table S3). Albeit being more abundant and overall shorter in length in RCC4 than in MDA-MB-231 cells, PADs showed a similar composition in the two cell lines. Namely, they were enriched in non-coding genes (Table S3) and in the constitutive heterochromatin mark H3K9me3 (Fig. 2A and B). In addition, by measuring the differential enrichment of repetitive DNA elements in PML ChIP-seq data with respect to input sequences in each cell line, we found that PADs are enriched in several classes of repetitive DNA elements such as LTRs, LINEs and DNA repeats in both MDA-MB-231 and RCC4 cells (Fig. 2C–Table S4). Importantly, although PADs in MDA-MB-231 and RCC4 cells generally share the same genomic features, the overlay of genomic sequences decorated by PML identified both common and unique PADs (Fig. 2D), indicating that PML associated to overlapping genomic regions as well as cell-specific chromatin across the genome of these cell lines.Fig. 2 PADs are structurally similar in MDA-MB-231 and RCC4 cells. (A) Genome browser view of PML and H3K9me3 Log2(ChIP/input) ratios (y axis range shown in brackets) and called PADs in MDA-MB-231 and RCC4 cells. For ChIP-seq of H3K9me3 data were downloaded from GSE226060 for MDA-MB-231 and from GSE143653 for A498 cells. The latter was used as representative of ccRCC cells. (B) Enrichment of H3K9me3 in PADs, random PADs (R-PADs) and inter PADs regions (I-PADs) identified in MDA-MB-231 and RCC4 cells; bar, median; whiskers, min-max; ****P < 10−4, unpaired t-tests with Welch's correction. (C) Enrichment plot of repetitive DNA elements in MDA-MB-231 (left) and RCC4 (right) cells. Each dot represents a repetitive DNA element, labels represent families of repetitive elements and colors represent repetitive DNA classes. Dashed lines mark -log10(AdjustedPvalues = 0.1). Log2FoldChanges were obtained comparing PML ChIP and input samples. (C) Venn diagram of overlapping PADs among MDA-MB-231 and RCC4 cells, Fisher's exact test, two-tailed. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 2

3.3 PADs are domains of cell-specific gene regulation by PML

In MDA-MB-231 cells we previously showed that PML regulates the epigenetic composition of PADs to facilitate cell-relevant transcriptional functions therein [7]. Specifically, we described PADs as discontinuous heterochromatic domains characterized at large by the constitutive heterochromatin mark H3K9me3 and containing small euchromatic domains of decreased PML association [7]. Interestingly, because we observed that PML regulates H3K9me3 deposition in PADs while at the same time promoting gene expression in the smaller euchromatic domains contained therein, we proposed that PML participates to the confinement of chromatin domains of transcriptional activity by regulating heterochromatin at surrounding regions [7]. To understand if PML regulates gene expression within PADs in a tissue-specific manner, we conducted a comparative analysis of the transcriptome regulated by PML in MDA-MB-231 and RCC4 cells and mapped PML-regulated genes to PADs (Table S6).

In line with our previous findings, gene set enrichment analysis revealed that in MDA-MB-231 cells the majority and most significantly regulated genes whose expression is promoted by PML pertain to cell motility processes, while negatively regulated genes are enriched in functional categories involved in protein transport (Fig. S3A). In RCC4 cells, PML had different regulatory effects on gene expression. Specifically, silencing of PML led to the downregulation of a large number of genes involved in cell cycle regulation and to the upregulation of genes clustering in immune-related defense responses (Fig. S3A). The overlap of PML regulated genes across these cell lines identified few common genes that clustered in functional categories distinct from the largest and most significant categories previously described (Figs. S3B and S3C). In line with these data, we defined cell-specific functions of PML towards promoting metastasis or cell cycle progression in TNBC versus ccRCC [12,13]. Therefore, these data indicate that PML exerts distinct oncogenic functions in different contexts by participating to the regulation of cell-specific transcriptional programs.

Next, we mapped the PML-regulated transcriptome to regions of PML DNA association. We confirmed that the large majority of PML-regulated genes contained within PADs in MDA-MB-231 cells are positively regulated by PML [7] and they are enriched in genes involved in cell motility process, such as TNC and ADAMTS9 (Fig. 3A and B, Table S6). Similarly, in RCC4 cells PML mostly promotes gene expression within PADs (Fig. 3C–Table S6), with genes belonging to cell cycle regulation gene sets, such as CDK1 (Fig. 3D). Of note, within all regulated genes contained within PADs, a minority are commonly regulated by PML in both cell lines while the majority shows cell-specific regulation (Fig. 3E). Interestingly, TFs motif analysis revealed an enrichment of KZFP (Kruppel-associated box (KRAB)-containing zinc finger proteins) predicted binding sites within commonly regulated genes, while cell-specific regulated genes are enriched in binding motifs of different transcription factors (Fig. 3F). These include known tissue-specific oncogenes like HIF2α (EPAS1), which is aberrantly expressed in ccRCC [14], and PROX1, which is overexpressed and promotes migration and invasion in TNBC [15]. In accordance with their function, predicted binding sites for these TFs were identified in genes regulated in RCC4 and MDA-MB-231 cells respectively (Fig. 3F). In sum, these data indicate that PADs are regions of cell-specific transcriptional activity promoted by PML.Fig. 3 PML mediates transcriptional regulation within PADs in MDA-MB-231 and RCC4 cells. (A) Enrichment plot of genes deregulated upon silencing of PML and falling into PADs in MDA-MB-231 cells. (B) Gene set enrichment analysis (GSEA) plot of a gene set containing a high number of genes downregulated upon PML silencing and mapping to PADs in MDA-MB-231 cells (vertical lines in the upper X-axis). The enrichment score (ES) on the y-axis reflects the degree of differential expression of the analyzed gene set (with negative ES reppresenting downregulated genes). Normalized enrichment score (NES) is calculated to determine the statistical significance of the ES. Associated P-values are represented. (C) Enrichment plot of genes deregulated upon silencing of PML and falling into PADs in RCC4 cells. (D) GSEA plot of a gene set containing a high number of genes downregulated upon PML silencing and mapping to PADs in RCC4 cells. Refer to panel B for data reppresentation. (E) Venn diagram of overlapping genes regulated by PML within PADs in RCC4 and MDA-MB-231 cells. (F) Top 5 TFs motifs enriched at 10 kb from the TSS of genes mapping to PADs and commonly or uniquely regulated by PML in RCC4 and MDA-MB-231. Data are represented as motif sequences, TF name and p-value.

Fig. 3

3.4 PML organizes common chromatin domains in cancer cell lines

To corroborate our data in another biological system, we used NB4 cells, a cellular model of acute promyelocytic leukemia (APL). APL is typified by the t(15; 17) chromosomal translocation that fuses the PML gene to the RARA gene, leading to the formation of the oncogenic PML-RARα fusion protein. Although the t(15; 17) translocation is monoallelic, PML-RARα interacts with and exerts dominant-negative effects on the remaining wild-type PML protein, leading to a genomic distribution of PML that is largely superimposable to that of PML-RARα [16]. Importantly, treatment of NB4 cells with all-trans retinoic acid (ATRA) induces degradation of PML-RARα and releases wild-type PML from it physical and functional sequestration [16]. Therefore, we took advantage of deposited PML ChIP-seq data from NB4 cells ± ATRA (GSE18886) to assess whether PML-RARα dislodged PML from PADs and ATRA treatment rescued PML genome association as defined thus far [16].

Analysis of PML chromatin association by EDD identified qualitatively and quantitatively different domains in NB4 cells treated or not with ATRA. Specifically, fewer and shorter PADs were identified in untreated NB4 cells, while ATRA treatment increased PADs numbers and PML genome coverage (Table S2 and Fig. 4A). Notably, measurement of repetitive elements in PML ChIP-seq versus input samples revealed that in untreated NB4 cells, PADs were depleted of repetitive elements, with the exception of satellite repeats at low significance, and ATRA treatment restored the enrichment of repetitive elements within PADs (Fig. 4B–Table S7). The classes of DNA repeats that became enriched in PML-associated chromatin upon ATRA treatment of NB4 cells are similar to repetitive elements found in PADs in MDA-MB-231 and RCC4 cells, namely LTRs, LINEs and DNA repeats (Fig. 4, Fig. 2B). Interestingly, 106 PADs were common between untreated and ATRA-treated NB4 cells, suggesting a potentially stronger association of PML to these regions. Alternatively, enrichment of satellite repeats, which represent highly heterochromatic regions [17], in PADs escaping the dominant negative effect of PML-RARα might be due to a high concentration of heterochromatic PML interactors that grant PML persistence. Common PADs contained a combined set of 459 genes, with a major contribution of genes belonging to the olfactory and KZFP gene clusters (Fig. 4C).Fig. 4 Analysis of PADs in NB4 cells untreated or treated with ATRA. (A) Genome browser view of PML Log2(ChIP/input) ratios (y axis range shown in brackets) and called PADs in untreated (-ATRA) and treated (+ATRA) NB4 cells. (B) Enrichment plot of repetitive DNA elements in untreated (-ATRA, left) and treated (+ATRA, right) NB4 cells. Each dot represents a repetitive DNA element, labels represent families of repetitive elements and colors represent repetitive DNA classes. Dashed lines mark -log10(AdjustedPvalues = 0.1). Log2FoldChanges were obtained comparing PML ChIP and input samples. Outliers including satellite repeats were removed from the plot to tidy the data. (C) Venn diagram of overlapping PADs and coding genes in untreated (-ATRA) and treated (+ATRA) NB4 cells (left), Fisher's exact test, two-tailed, and KEGG enriched pathways of overlapping genes (right). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 4

Taken together, these findings provide further evidence to the dominant negative effect of PML-RARα over the wild-type PML protein in APL [16] by showing a delocalization of PML from genomic domains enriched in repetitive elements, which are restored upon ATRA treatment.

The analysis of PML chromatin association in NB4 cells gave us the opportunity to further compare the genetic identity of PADs across cell lines by overlapping PML-associated DNA in MDA-MB-231, RCC4 and NB4 cells ± ATRA. When using untreated NB4 cells, a total of 72 PADs were common to the 3 cell lines (Fig. S4A). We defined these as common PADs (cPADs) and confirmed that they are enriched in olfactory and KZFP gene clusters (Fig. S4B). These data suggest that PML associates constitutively to these regions, irrespective of cell identity and even in conditions of PML functional inactivation. Accordingly, PML association to similar gene clusters was also reported in MEFs [6], indicating that this genomic association is also species independent.

When comparing MDA-MB-231, RCC4, and NB4 cells + ATRA, a larger number of cPADs were identified, containing a greater number of genes (Fig. 5A, Table S8). Gene set enrichment analysis showed enrichment of several duplicated gene clusters, that is families of genes with conserved functional activity and proximal genomic localization, such as the ATP-binding cassette (ABC) transporters and amylase genes in addition to the olfactory and KZFP gene clusters (Fig. 5B). Consistently, analysis of cPADs gene distribution across the human genome showed that the majority of these genes clustered within spatially restricted regions (Fig. 5C–Table S9), in agreement with the overrepresentation of duplicated gene clusters. Finally, TF motif calling within cPADs revealed an enrichment for KZFP (Fig. 5D), in line with the presence of repetitive DNA in these regions [17]. Of note, the vast majority of genes contained within cPADs are not transcriptionally regulated by PML in MDA-MB-231 or RCC4 cells (PML regulates 6 % and 9 % of cPADs genes in MDA-MB-231 and RCC4 cells respectively), in line with their cell type-restricted expression [18].Fig. 5 Common PADs are enriched in duplicated gene clusters. (A) Venn diagram of overlapping PADs and coding genes in RCC4, MDA-MB-231 and ATRA treated (+ATRA) NB4 cells. Chi-squared test, p-value<0,01. (B) KEGG enriched pathways of coding genes mapping to common PADs. (C) Karyotype plot of coding genes, represented by red dots, mapping to common PADs (cPADs). Blue-filled dots indicate genes that are statistically enriched in clusters. (D) Top 5 TF motifs enriched in cPADs. Data are represented as motif sequence, TF name and p-value. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 5

In summary, these data unveil a recurrent pattern of PML chromatin association that is centered on repetitive elements and duplicated gene clusters.

4 Discussion

In this study we applied a computational protocol to identify PML-associated chromatin in ChIP-seq experiments and compared PML chromatin occupancy across different cancer cell lines. In so doing, we identified common and cell specific PML chromatin association. To our knowledge, this is the first study that unveils a modality of PML chromatin association that is similar across different cancer cell lines. We reveal that PML associates to repetitive DNA elements and duplicated gene clusters that reside in constitutive heterochromatin outside embryonal cells or specific cell lineages [18]. Therefore, our findings provide support to a growing body of evidence that implicates PML as an important regulator of heterochromatin organization [2].

PML reportedly exerts complex molecular functions of transcriptional regulation, which have been described as direct vs indirect and originating from the PML-NBs of from nucleoplasmic PML moieties (Fig. 6). To obtain a rigorous definition of PML chromatin association regardless of its localization to the PML-NBs, we first took a comparative computational approach to map PML chromatin association upon ChIP-seq in a cell line of TNBC. By comparing 4 peak calling algorithms based on different computational models we could measure PML DNA association with statistical confidence only when using EDD, an algorithm that enables detection of chromatin association by proteins that are widely distributed and with a low level of enrichment on DNA, which was previously used to define large heterochromatic PML-associated domain or PADs [8]. This type of chromatin is fundamentally different from the DNA that was identified by adjacency to the PML-NBs, which is enriched in euchromatic gene regulatory elements of narrow PML association [2]. A logical implication of our data is that chromatin immunoprecipitation protocols enrich nucleoplasmic PML moieties, a tenet that is confirmed by lack of association of PADs with the PML-NBs [6,7]. Therefore, with our comparative analysis we provide a methodological approach to analyze PML ChIP-seq experiments.Fig. 6 Molecular functions of PML in transcriptional regulation. Nuclear PML is distributed in nucleoplasmic moieties and large aggregates known as PML-nuclear bodies (PML-NBs). Vast literature has demonstrated that transcription factors (purple) and transcriptional regulators like epigenetic factors (green) transit through the PML-NBs to be regulated. Emerging literature is suggesting that nucleoplasmic PML associates with and regulates large heterochromatic regions (gray) that are enriched in repetitive elements (yellow) and duplicated gene clusters (dark blue). Biorender.com was used to create the figure. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 6

Having defined a standardized way for analyzing PML ChIP-seq, we asked whether PML chromatin association in PADs is similar in other cancer cell lines besides TNBC, and we selected two additional cell lines where PML exerts important, albeit opposite functions: an oncogenic role in ccRCC cells [12] and a tumor suppressive function inhibited by the dominant negative oncogene PML-RARα in APL cells [16]. In these contexts, we found that a sizable number of PADs are similar across cell lines, along with the presence of cell type-specific PADs, and that PADs share the same genomic features, that is enrichment in non-coding genes and repetitive DNA elements mostly of the LINE and LTR families. However, PADs have on average different lengths in the different cell lines, which may be explained by genetic/epigenetic differences but also by a differential enrichment of PML on chromatin. In this respect, we have recently demonstrated that the relative concentration of PML at PML-NBs vs nucleoplasm is higher in RCC4 than in MDA-MB-231 cells [12]. Because we previously showed that PADs are at large formed by nucleoplasmic PML [7], it could be speculated that PADs identified in RCC4 cells are shorter because of lower amounts of nucleoplasmic PML.

Of note, PML association to PADs is reduced in cells that express the PML dominant negative protein PML-RARα and becomes restored upon reestablishment of wild-type PML function. These data indicate that repetitive elements are constitutive components of PADs and suggest that the association of PML to DNA is intricately liked to the distribution of DNA repeats. Of note, the association of PML at common PADs may appear at odds with the differential oncogenic vs tumor suppressive function of PML in the cell types that we have analyzed. However, PML is a complex protein with distinct biochemical (i.e. via the PML-NBs or nucleoplasmic PML) and biological functions, and being PML association to PADs one amongst many functions of PML, it is plausible that the differential activities of PML in these cellular contexts are not mediated by common PADs.

Nonetheless, an interesting observation stems from the subset of PADs that are common across the cell lines that we analyzed, which contain functionally related genes. These include the olfactory, KZFPs, GABA receptor and ABC transporter gene families, all of which are characterized by their organization in gene clusters, heterochromatinization and cell type restricted expression [18]. Mechanistically, the establishment of heterochromatin at functionally related gene clusters and repetitive elements is driven by KZFPs, which recruit on DNA a repressive complex orchestrated by KRAB-associated protein 1 (KAP1/TRIM28) and containing the histone methyltransferase SETDB1 and additional chromatin regulators like heterochromatin protein HP1α [19]. Motif calling analysis showed that common PADs are enriched in KZFPs binding motifs. Because it was recently shown that PML interacts with KAP1 and SETDB1 and promotes their DNA association [7,20], taken together these data suggest that PML may functionally regulate the activity of a multiprotein repressor complex that promotes epigenetic silencing of repetitive and cluster DNA regions. Albeit PML association to PADs reportedly occurs at least in part outside the PML-NBs, PML association to KAP1 occurs at the PML-NBs, resulting in its SUMOylation and activation [20]. Along similar lines, KZFPs-KAP1 nuclear foci were described as adjacent to the PML-NBs in embryonic cells [21]. In this scenario, it is tempting to speculate that PML partakes to the assembly of this repressive complex at the PML-NBs, a complex that is then recruited, along with PML moieties, to repetitive DNA elements. Further experiments will be necessary to test this hypothesis.

Interestingly, in addition to LINE and LTR elements, PADs are enriched in the class of DNA repeats (such as PiggyBac and Mariner transposons), which are amongst the oldest repetitive elements. Transposons belonging to this class reproduce via horizontal transfer, a known virulence mechanism [22] and have been latent for the past 40 Mya [23]. Based on their nature, it is thus tempting to speculate that the association of PML to their genomic regions may have participated to their suppression during evolution. This hypothesis is in line with an interesting evolutionary link that was recently established between the molecular functions of PML and retroelements. It was suggested that the PML gene first emerged during evolution to suppress L1 retrotransposition by providing a cytoplasmic exonuclease activity [24]. Similarly, KZFPs are evolutionarily linked to repetitive elements of ancient viral origin, like transposable elements, as a means to restrict their uncontrolled activation and ensuing chromosomal instability [25]. Thus, nuclear cooptation of PML functions may have occurred along evolution to widen its co-repressive functions over impending genome occupation by repetitive elements and may position PML at an evolutionary crossroad in the epigenetic control of chromatin instability.

Limitations of the study

This study has some limitations. First, we have compared only 4 peak calling algorithms and cannot exclude that there are additional computational methods that work as well or better than EDD to define PML-associated chromatin. Second, we have analyzed only 3 cell lines of different tumor origin and did not characterize primary human cells. Third, we did not provide evidence that PML promotes deposition of repressive epigenetic modifications at all PADs in all cell lines. Fourth, we do not provide functional evidence of the role of PML in restricting gene expression or retrotransposition from repetitive DNA elements. Nonetheless, our data provide a starting point for several lines of investigation. These include delving into the evolutionary co-option of PML-like genes along with KZFPs towards the epigenetic repression of transposable elements and the potential transcriptional regulation of repetitive DNA elements by PML.

Resource availability

Lead contact

Requests for further information and reagents should be directed to and will be fulfilled by the Lead Contact, Doctor Rosa Bernardi (bernardi.rosa@hsr.it).

Data availability

Original data associated with this study have been deposited into a publicly available repository. Specifically, PML ChIP-seq data performed in RCC4 cells were deposited in GEO under the accession number GSE261929. All other sequencing data were obtained from the publicly available repository GEO with the following accession numbers: PML ChIP-seq data from MDA-MB-231 cells: GSE226060; PML ChIP-seq data from NB4 cells: GSE18886; input samples from RCC4 cells: GSE120887; H3K9me3 ChIP-seq data from A498 cells: GSE143653; H3K9me3 ChIP-seq in MDA-MB-231 cells: GSE226060; RNA-seq data from MDA-MB-231 cells: GSE226111; RNA-seq data from RCC4 cells: GSE246846.

CRediT authorship contribution statement

Cristina Fracassi: Writing – original draft, Validation, Methodology, Investigation, Formal analysis, Data curation. Matilde Simoni: Investigation. Martina Uggè: Investigation. Marco J. Morelli: Methodology, Data curation, Conceptualization. Rosa Bernardi: Writing – original draft, Supervision, Funding acquisition, Conceptualization.

Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:Rosa Bernardi reports financial support was provided by Italian Association for Cancer Research. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A Supplementary data

The following are the Supplementary data to this article.

Acknowledgements

The study was funded by 10.13039/501100005010 Italian Association for Cancer Research (AIRC) with an Investigator Grant [20170 to R.B.].

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.heliyon.2024.e36499.
==== Refs
References

1 Abou-Ghali M. Lallemand-Breitenbach V. PML Nuclear bodies: the cancer connection and beyond Nucleus 15 2024 1 13 10.1080/19491034.2024.2321265
2 Corpet A. Kleijwegt C. Roubille S. Juillard F. Jacquet K. Texier P. Lomonte P. PML nuclear bodies and chromatin dynamics: catch me if you can Nucleic Acids Res. 48 2020 11890 11912 10.1093/nar/gkaa828 33068409
3 Ching R.W. Ahmed K. Boutros P.C. Penn L.Z. Bazett-Jones D.P. Identifying gene locus associations with promyelocytic leukemia nuclear bodies using immuno-TRAP J. Cell Biol. 201 2013 325 335 10.1083/jcb.201211097 23589495
4 Wang M. Wang L. Qian M. Tang X. Liu Z. Lai Y. Ao Y. Huang Y. Meng Y. Shi L. Peng L. Cao X. Wang Z. Qin B. Liu B. PML2‐mediated thread‐like nuclear bodies mark late senescence in Hutchinson–Gilford progeria syndrome Aging Cell 19 2020 1 14 10.1111/acel.13147
5 Kurihara M. Kato K. Sanbo C. Shigenobu S. Ohkawa Y. Fuchigami T. Miyanari Y. Genomic profiling by ALaP-seq reveals transcriptional regulation by PML bodies through DNMT3A exclusion Mol. Cell 78 2020 493 505.e8 10.1016/j.molcel.2020.04.004 32353257
6 Delbarre E. Ivanauskiene K. Spirkoski J. Shah A. Vekterud K. Moskaug J.Ø. Bøe S.O. Wong L.H. Küntziger T. Collas P. PML protein organizes heterochromatin domains where it regulates histone H3.3 deposition by ATRX/DAXX Genome Res. 27 2017 913 921 10.1101/gr.215830.116 28341773
7 C. Fracassi, M. Ugge’, M. Abdelhalim, E. Zapparoli, M. Simoni, D. Magliulo, D. Mazza, D. Lazarevic, M.J. Morelli, P. Collas, R. Bernardi, PML modulates epigenetic composition of chromatin to regulate expression of pro-metastatic genes in triple-negative breast cancer, Nucleic Acids Res. 51 (2023) 11024–11039. 10.1093/nar/gkad819.
8 Lund E. Oldenburg A.R. Collas P. Enriched domain detector: a program for detection of wide genomic enrichment domains robust against local variations Nucleic Acids Res. 42 2014 10.1093/nar/gku324 e92–e92
9 Xu S. Grullon S. Ge K. Peng W. Spatial clustering for identification of ChIP-enriched regions (SICER) to map regions of histone methylation patterns in embryonic stem cells Kidder B.L. Methods Mol Biol 2014 Springer New York New York, NY 97 111 10.1007/978-1-4939-0512-6_5
10 Zhang Y. Liu T. Meyer C.A. Eeckhoute J. Johnson D.S. Bernstein B.E. Nusbaum C. Myers R.M. Brown M. Li W. Liu X.S. Model-based analysis of ChIP-seq (MACS) Genome Biol. 9 2008 R137 10.1186/gb-2008-9-9-r137 18798982
11 Feng X. Grossman R. Stein L. PeakRanger: a cloud-enabled peak caller for ChIP-seq data BMC Bioinf. 12 2011 139 10.1186/1471-2105-12-139
12 Simoni M. Menegazzi C. Fracassi C. Biffi C.C. Genova F. Tenace N.P. Lucianò R. Raimondi A. Tacchetti C. Brugarolas J. Mazza D. Bernardi R. PML restrains p53 activity and cellular senescence in clear cell renal cell carcinoma EMBO Mol. Med. 2024 10.1038/s44321-024-00077-3
13 Ponente M. Campanini L. Cuttano R. Piunti A. Delledonne G.A. Coltella N. Valsecchi R. Villa A. Cavallaro U. Pattini L. Doglioni C. Bernardi R. PML promotes metastasis of triple-negative breast cancer through transcriptional regulation of HIF1A target genes JCI Insight 2 2017 1 15 10.1172/jci.insight.87380
14 Hoefflin R. Harlander S. Schäfer S. Metzger P. Kuo F. Schönenberger D. Adlesic M. Peighambari A. Seidel P. Chen C. Consenza-Contreras M. Jud A. Lahrmann B. Grabe N. Heide D. Uhl F.M. Chan T.A. Duyster J. Zeiser R. Schell C. Heikenwalder M. Schilling O. Hakimi A.A. Boerries M. Frew I.J. HIF-1α and HIF-2α differently regulate tumour development and inflammation of clear cell renal cell carcinoma in mice Nat. Commun. 11 2020 4111 10.1038/s41467-020-17873-3 32807776
15 Zhu L. Tian Q. Gao H. Wu K. Wang B. Ge G. Jiang S. Wang K. Zhou C. He J. Liu P. Ren Y. Wang B. PROX1 promotes breast cancer invasion and metastasis through WNT/β-catenin pathway via interacting with hnRNPK Int. J. Biol. Sci. 18 2022 2032 2046 10.7150/ijbs.68960 35342346
16 Martens J.H.A. Brinkman A.B. Simmer F. Francoijs K.-J. Nebbioso A. Ferrara F. Altucci L. Stunnenberg H.G. PML-RARα/RXR alters the epigenetic landscape in acute promyelocytic leukemia Cancer Cell 17 2010 173 185 10.1016/j.ccr.2009.12.042 20159609
17 Haws S.A. Simandi Z. Barnett R.J. Phillips-Cremins J.E. 3D genome, on repeat: higher-order folding principles of the heterochromatinized repetitive genome Cell 185 2022 2690 2707 10.1016/j.cell.2022.06.052 35868274
18 Lu J.Y. Shao W. Chang L. Yin Y. Li T. Zhang H. Hong Y. Percharde M. Guo L. Wu Z. Liu L. Liu W. Yan P. Ramalho-Santos M. Sun Y. Shen X. Genomic repeats categorize genes with distinct functions for orchestrated regulation Cell Rep. 30 2020 3296 3311.e5 10.1016/j.celrep.2020.02.048 32160538
19 Bersaglieri C. Kresoja-Rakic J. Gupta S. Bär D. Kuzyakiv R. Panatta M. Santoro R. Genome-wide maps of nucleolus interactions reveal distinct layers of repressive chromatin domains Nat. Commun. 13 2022 1483 10.1038/s41467-022-29146-2 35304483
20 Tessier S. Ferhi O. Geoffroy M.-C. González-Prieto R. Canat A. Quentin S. Pla M. Niwa-Kawakita M. Bercier P. Rérolle D. Tirard M. Therizols P. Fabre E. Vertegaal A.C.O. de Thé H. Lallemand-Breitenbach V. Exploration of nuclear body-enhanced sumoylation reveals that PML represses 2-cell features of embryonic stem cells Nat. Commun. 13 2022 5726 10.1038/s41467-022-33147-6 36175410
21 Briers S. Crawford C. Bickmore W.A. Sutherland H.G. KRAB zinc-finger proteins localise to novel KAP1-containing foci that are adjacent to PML nuclear bodies J. Cell Sci. 122 2009 937 946 10.1242/jcs.034793 19258395
22 Emamalipour M. Seidi K. Zununi Vahed S. Jahanban-Esfahlan A. Jaymand M. Majdi H. Amoozgar Z. Chitkushev L.T. Javaheri T. Jahanban-Esfahlan R. Zare P. Horizontal gene transfer: from evolutionary flexibility to disease progression Front. Cell Dev. Biol. 8 2020 10.3389/fcell.2020.00229
23 Pace J.K. Feschotte C. The evolutionary history of human DNA transposons: evidence for intense activity in the primate lineage Genome Res. 17 2007 422 432 10.1101/gr.5826307 17339369
24 Mathavarajah S. Vergunst K.L. Habib E.B. Williams S.K. He R. Maliougina M. Park M. Salsman J. Roy S. Braasch I. Roger A.J. Langelaan D.N. Dellaire G. PML and PML-like exonucleases restrict retrotransposons in jawed vertebrates Nucleic Acids Res. 51 2023 3185 3204 10.1093/nar/gkad152 36912092
25 Helleboid P. Heusel M. Duc J. Piot C. Thorball C.W. Coluccio A. Pontis J. Imbeault M. Turelli P. Aebersold R. Trono D. The interactome of <scp>KRAB</scp> zinc finger proteins reveals the evolutionary history of their functional diversification EMBO J. 38 2019 1 16 10.15252/embj.2018101220
26 Dull T. Zufferey R. Kelly M. Mandel R.J. Nguyen M. Trono D. Naldini L. A third-generation lentivirus vector with a conditional packaging system J. Virol. 72 1998 8463 8471 10.1128/JVI.72.11.8463-8471.1998 9765382
27 Cabianca D.S. Casa V. Bodega B. Xynos A. Ginelli E. Tanaka Y. Gabellini D. A long ncRNA links copy number variation to a polycomb/trithorax epigenetic switch in FSHD muscular dystrophy Cell 149 2012 819 831 10.1016/j.cell.2012.03.035 22541069
28 Li H. Durbin R. Fast and accurate long-read alignment with Burrows–Wheeler transform Bioinformatics 26 2010 589 595 10.1093/bioinformatics/btp698 20080505
29 Amemiya H.M. Kundaje A. Boyle A.P. The ENCODE blacklist: identification of problematic regions of the genome Sci. Rep. 9 2019 9354 10.1038/s41598-019-45839-z 31249361
30 Criscione S.W. Zhang Y. Thompson W. Sedivy J.M. Neretti N. Transcriptional landscape of repetitive elements in normal and cancer human cells BMC Genom. 15 2014 583 10.1186/1471-2164-15-583
31 Langmead B. Aligning short sequencing reads with Bowtie Curr. Protoc. Bioinforma. 32 2010 1 24 10.1002/0471250953.bi1107s32
32 Love M.I. Huber W. Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 Genome Biol. 15 2014 550 10.1186/s13059-014-0550-8 25516281
33 Ginestet C. ggplot2: elegant graphics for data analysis J. R. Stat. Soc. Ser. A (Statistics Soc. 174 2011 245 246 10.1111/j.1467-985X.2010.00676_9.x
34 Bolger A.M. Lohse M. Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data Bioinformatics 30 2014 2114 2120 10.1093/bioinformatics/btu170 24695404
35 Dobin A. Davis C.A. Schlesinger F. Drenkow J. Zaleski C. Jha S. Batut P. Chaisson M. Gingeras T.R. STAR: ultrafast universal RNA-seq aligner Bioinformatics 29 2013 15 21 10.1093/bioinformatics/bts635 23104886
36 Liao Y. Smyth G.K. Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features Bioinformatics 30 2014 923 930 10.1093/bioinformatics/btt656 24227677
37 Quinlan A.R. Hall I.M. BEDTools: a flexible suite of utilities for comparing genomic features Bioinformatics 26 2010 841 842 10.1093/bioinformatics/btq033 20110278
38 Neph S. Kuehn M.S. Reynolds A.P. Haugen E. Thurman R.E. Johnson A.K. Rynes E. Maurano M.T. Vierstra J. Thomas S. Sandstrom R. Humbert R. Stamatoyannopoulos J.A. BEDOPS: high-performance genomic feature operations Bioinformatics 28 2012 1919 1920 10.1093/bioinformatics/bts277 22576172
39 Khan A. Mathelier A. Intervene: a tool for intersection and visualization of multiple gene or genomic region sets BMC Bioinf. 18 2017 287 10.1186/s12859-017-1708-7
40 Ramírez F. Dündar F. Diehl S. Grüning B.A. Manke T. deepTools: a flexible platform for exploring deep-sequencing data Nucleic Acids Res. 42 2014 W187 W191 10.1093/nar/gku365 24799436
41 Robinson J.T. Thorvaldsdóttir H. Winckler W. Guttman M. Lander E.S. Getz G. Mesirov J.P. Integrative genomics viewer Nat. Biotechnol. 29 2011 24 26 10.1038/nbt.1754 21221095
42 Heinz S. Benner C. Spann N. Bertolino E. Lin Y.C. Laslo P. Cheng J.X. Murre C. Singh H. Glass C.K. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities Mol. Cell 38 2010 576 589 10.1016/j.molcel.2010.05.004 20513432
43 Korotkevich G. Sukhov V. Budin N. Atryomov M.N. Sergushichev A. Fast Gene Set Enrichment Analysis 2021 BioRxiv 1 29 10.1101/060012 bioRxiv
44 Ge S.X. Jung D. Yao R. ShinyGO: a graphical gene-set enrichment tool for animals and plants Bioinformatics 36 2020 2628 2629 10.1093/bioinformatics/btz931 31882993
45 Gel B. Serra E. karyoploteR: an R/Bioconductor package to plot customizable genomes displaying arbitrary data Bioinformatics 33 2017 3088 3090 10.1093/bioinformatics/btx346 28575171
46 Pazos Obregón F. Soto P. Lavín J.L. Cortázar A.R. Barrio R. Aransay A.M. Cantera R. Cluster Locator, online analysis and visualization of gene clustering Bioinformatics 34 2018 3377 3379 10.1093/bioinformatics/bty336 29701747
