HiCap depicts interactome dynamics of blood stem-like cells upon drug exposure

We processed 0.8 and 1.1 billion reads mapped on the transcriptome and interactome of 18 samples corresponding to normal and drug-treated counterparts of the three cell lines: CMK, MOLM-1, and K-562, in two replicates. These cell lines are a good and convenient surrogate model readily available in many research laboratories (Skopek et al, 2023). We first investigated gene expression changes upon drug exposure. The technical replicates of the experiments correlated well in all cases (coefficient of correlation >0.93 across all conditions) (Fig S1A–F). However, RNA profiles of samples were not always able to separate treated cells from their normal counterparts (Fig S2). Nevertheless, we detected differentially expressed (DE) genes (FDR < 0.1), absolute log-fold change (abs[logFC] > 1.2) belonging to relevant processes for drug exposure such as regulation of hemopoiesis (GO:1903706, FDR = 1 × 10⁻⁸) and response to toxic substances (GO:0009636, FDR = 5 × 10⁻⁴) (Fig S3A and B, Table S1), confirming the cells’ response to the treatment.

Next, we analyzed HiCap datasets to find genomic interactions of targeted promoters and selected variants (see the Materials and Methods section). We detected 114,365 distal regions (i.e., untargeted regions that are located distal to the promoters) interacting with 7,254 promoters (supporting pairs >4 and Bonferroni-adjusted P-value < 0.01) across all cell lines and treatments (Table S2). Promoter-interacting regions (PIRs) covered around 2.8% (84.4 Mb) of the genome. We also detected 18,024 interactions between promoters and 19,204 interactions between targeted variants and distal regions. Because K-562 is a tier 1 ENCODE cell line, a vast array of public ATAC-seq and ChIP-seq datasets is available (Table S3). We used these datasets to assess the enhancer potential of the PIRs. In K-562 cells, 51.2% (22,719/44,366) of the PIRs overlapped with at least one enhancer mark (6.52-fold enrichment) (Fig 1A). We found that 41% (18,190/44,366) of the PIRs overlapped only a TF-binding site, highlighting the importance of TF datasets during the evaluation of promoter–enhancer interactions. All super-enhancers (SEs) annotated in K-562 cells (Jiang et al, 2019) also overlapped with at least one PIR.

We further tested the potential of detected PIRs to modulate the activity of targeted genes using an inducible dual-effector epigenetic interference system enCRISPRi-LK (Li et al, 2020). We successfully reduced the gene expression in two out of three cases by repressing their enhancers (Fig 1B and C). The first gene, TGFBR1 (Transforming Growth Factor β Receptor 1), is a part of the transforming growth factor β signaling pathway known for supporting the maintenance of the self-renewal capacity of hematopoietic stem cells (Blank & Karlsson, 2015). It is regulated by an enhancer located 272 kb away in the chr9:101,595,333–101,595,751 region (Fig S4A). The second gene, ANAPC4 (Anaphase Promoting Complex Subunit 4), is part of a highly conserved protein complex that controls the amounts of the cyclins and other cell cycle regulators to ensure proper cell cycle transitions (Wäsch et al, 2010). Enhancer of ANAPC4 is located 326 kb away from it in the chr4:25,051,758–25,052,804 region (Fig S4B). None of these regions interact with NFE2 (Nuclear Factor Erythroid 2) promoter used as a negative control.

We then investigated different aspects of interaction data by looking at interaction statistics, functional enrichment of interacting regions, and expression profiles of interacting promoters. The average interaction distance across all the samples is 156 kb (Fig 2A). Our interactome maps showed close to single-enhancer resolution (average PIR length of 860 bp) and therefore facilitated the individual discovery of shared, and cell type-specific interactions (Fig 2B). Interactome profiles were more effective in separating untreated cells from their treated counterparts compared with that of the transcriptome data (Fig S2). Interestingly, 1,302 (43.9%) of targeted variants showed at least one interaction in one cell type, and interacting variants had significantly lower minor allele frequencies compared with those targeted variants without an interaction (Fig 2C, minor allele frequencies = 0.16 versus 0.26, FDR = 7.58 × 10⁻¹⁷¹).

The ChIP-seq profiles of 218 TFs from K-562 cells (Table S3) were used to investigate the TF distribution across PIRs. The number of TF-binding sites per PIR followed an inverse power–law distribution; that is, most PIRs contained few TFs, and few PIRs contained many TFs (Fig 2D). Using the k-means clustering algorithm, we observed four main TF clusters that co-occur, whereas generally the binding patterns of most TFs were not correlated (Fig S5) suggesting a diversity of enhancers based on their TF-binding patterns. One of the clusters contained 49 TFs that are highly enriched for SWI/SNF superfamily (FDR = 6.3 × 10⁻¹⁶), NuRD complex (FDR = 3.8 × 10⁻¹³), and leukemogenesis (FDR = 5.8 × 10⁻¹¹). That is in line with previous studies showing a high mutation rate of chromatin remodeling complexes in cancer (Bracken et al, 2019). We then separated PIRs connected to either expressed (top 75% quantile, 16,092 interactions) or not-expressed (bottom 25% quantile, 16,241 interactions) genes to see if there are any differential TF binding between the two groups. There were 52 and 24 TFs binding to PIRs connecting to high or low-expressed genes, respectively, and 23 TF were shared (Table S4). Transcriptional repressor DEAF1 (Michelson et al, 1999) was only found in PIRs connected to low-expressed genes. Moreover, 29 TFs were only found in PIRs connected to expressed genes (Fig S6A and B). Among them, members of the histone deacetylase family (HDAC1, HDAC2, EP400, NCOR1, and TBL1XR1) were present, confirming the role of histone deacetylation in gene activation (Jepsen et al, 2000; Bhaskara et al, 2008; LaBonte et al, 2009; Wang et al, 2009; Seto & Yoshida, 2014). Lastly, we looked at the TF-binding profiles of PIRs as a function of their overlap with SEs. Members of activator complexes such as EP300, EP400, TBL1XR1, NR2C2, and SOX6 were only found in PIRs overlapping with SEs (Fig S6C). However, CTCF and RAD21 (subunit of the cohesin complex [Hauf et al, 2001]) were found only in PIRs not overlapping with SEs (Fig S6D), which supports their role in regulating loop stability (Hansen et al, 2017).

Finally, we focused on the targeted genes across all cell lines and treatments. There were 8,251, 6,582, and 8,142 interacting genes in CMK, MOLM-1, and K-562 cells, respectively. We discovered cell type-specific and treatment-specific interacting genes (Figs 2E and S7A and B). In all cells, they were more likely to be expressed (Fig 2F). However, they were not more likely to be differentially expressed, supporting the different dynamics of the transcriptional regulation by enhancers and promoters (Larsson et al, 2019).

Network analysis of interactome facilitates the discovery of differentially interacting genes

To track the dynamics of promoter–enhancer interactions, we constructed networks using promoters and PIRs as network nodes and connected them with an edge in case of interaction (see the Materials and Methods section). We generated a separate network for each replicate of cell type and state. We then calculated two parameters, namely, Overlap Coefficient (OCE) and Jaccard Similarity Index (JI), to quantify the connectivity differences for each node across different states for each cell type. If the OCE/JI is equal to one, it means that the gene did not change its interactions. In contrast, if the OCE/JI equals zero, then it means that the gene either changed all interactions, gained all new interaction(s) or lost all its interaction(s) upon drug exposure. As an example, EEFG1 (Eukaryotic Translation Elongation Factor 1 Gamma) is a housekeeping gene involved in the translation mechanisms (Kumabe et al, 1992), and its interaction profile did not change significantly in CMK cells upon carboplatin induction (Fig 3A). Meanwhile, EYA3 (EYA Transcriptional Coactivator And Phosphatase 3) plays a particular function as a distinguishing mark between apoptotic and repair responses to genotoxic stress (Cook et al, 2009; Krishnan et al, 2009). It entirely changed the interaction profile under the same conditions (Fig 3B). We named the promoter nodes that changed their connectivity as differentially interacting (DI) genes (see the Materials and Methods section). This approach revealed hundreds of genes with interaction changes upon treatment despite their relatively even steady-state expression levels (Table S5).

In CMK cells, there were 156/45 DE genes and 728/696 DI genes upon carboplatin and gemcitabine treatments, respectively. Generally, in CMK cells, only DI genes showed enrichments for hematological cell count traits. Therefore, we assessed the biological relevance of DI genes by comparing their enrichment for human phenotype terms with that of non-DI genes, that is, genes that did not change their interactions upon treatment. Indeed, DI genes in CMK cells were more enriched for relevant phenotypes (Fig 3C). For example, ETV6 (ETS Variant Transcription Factor 6) is a transcriptional repressor implicated in dominantly inherited thrombocytopenia (Hock & Shimamura, 2017). Despite no significant changes in steady-state expression levels, it underwent considerable interactome changes upon both carboplatin and gemcitabine exposure (Fig 3D).

In MOLM-1 cells, we detected 11/252 DE genes and 464/625 DI genes upon exposure to carboplatin/gemcitabine, respectively. They were enriched for relevant traits except for DE genes upon exposure to carboplatin (Fig S8A and B). In particular, MPL (Proto-Oncogene, Thrombopoietin Receptor) is essential for the proliferation of megakaryocytes and platelet differentiation (Ng et al, 2014). It was not differentially expressed, but it changed its interaction landscape significantly upon exposure to both drugs (Fig S8C).

In K-562 cells, there were 555/44 DE genes and 921/932 DI genes upon treatment with carboplatin/gemcitabine, respectively. DE genes, upon carboplatin treatment, were enriched for multiple processes related to cell cycle regulation and hematological processes, whereas DI genes showed limited enrichments (Fig S9A). Conversely, very few genes were differentially expressed upon gemcitabine treatment, but DI genes were enriched for multiple processes related to lymphocyte and myeloid cell processes (Fig S9B). ERG (ETS Transcription Factor ERG), an essential gene for hematopoiesis (Knudsen et al, 2015), changed both its steady-state expression and interaction status in both treatments (Fig S9C). Concordantly, its haploinsufficiency impairs the self-renewal of hematopoietic stem cells under myelotoxic stress (Ng et al, 2011). The summary of DE and DI genes across cell types and treatments is presented in Table 1.

Genes connected to candidate enhancer variants tend to differentially interact upon drug exposure

In our previous study, we characterized 8,072,672 single-nucleotide variants (SNVs) from 96 NSCLC patients with differing chemotherapy-induced myelosuppression levels. Only a few protein-coding mutations were associated with the toxicity response of the patients (Björn et al, 2020a). We hypothesized that some SNVs might modulate the activity of enhancers involved in myelosuppression or entirely disrupt them. We overlapped SNVs with 95,567 PIRs derived from HiCap experiments on model cell lines used here to test our hypothesis. Half of the PIRs (52% or 49,802 PIRs) contained at least one variant (50,119 SNVs), and 94.2% of those contained only a single variant.

To filter SNVs relevant to higher or lower risk of drug toxicity, we grouped variants based on their allele count (AC) and allele frequency (AF) differences between the two patient cohorts: high toxicity (HT) and low toxicity (LT) (Fig S10). Accordingly, there were 6,583 variants whose AC and AF differed at least by 25% (see the Materials and Methods section). We also required these variants to have a strong effect size (ES) on TF binding affinity. Using the motifbreakR package, we calculated the ES of the selected variants across a collection of PWMs. Moreover, we filtered only cases where both the interacting gene and TF were expressed (mean [transcripts per million (TPM) across the cell line under all treatments] > 0.2). That prioritized the top 2,720 variants further considered as candidate variants.

Afterward, we asked if promoters connected to PIRs carrying these candidate variants are likelier to change their interaction profile upon drug treatment than those connected to PIRs carrying other (rest) variants. Indeed, promoters interacting with PIRs carrying candidate variants were more likely to change their interactions upon exposure to carboplatin (Fig 4A–C) and gemcitabine (Fig 4D–F). The mean fraction of DI genes carrying variants in PIRs varied from 5.0–7.2% for candidate variants to 1.6–2.7% for the rest cases. The highest P-value between the two groups was 1.7 × 10⁻⁵ in MOLM-1 cells in normal versus carboplatin-treated cases (Fig 4B), proving that interactions bearing candidate variants are significantly enriched for DI genes. As a result, we report 196 candidate DI genes connected to the candidate variants, which may be involved in the genetic mechanisms of chemotherapy-induced myelosuppression (Table S6).

We further investigated the relevance of the candidate gene set for the myelosuppression. We isolated genes associated (through mutations overlapping genes only) with the counts of lymphocytes, erythrocytes, platelets, and neutrophils from all associations NHGRI-EBI GWAS Catalogue v.1.0.2 (Sollis et al, 2023). Candidate genes showed a significant (P-value = 1.4 × 10⁻³) enrichment for blood cells count traits.

Finally, we investigated which TFs are disrupted by candidate variants in the enhancers of candidate genes. They showed significant enrichment for abnormal blood cell morphology/development (Bonferroni-adjusted P-value = 6.881 × 10⁻²⁹) and other related mouse phenotype terms (Table S6). Moreover, they broadly overlapped with TFs previously associated with NSCLC and different stages of chemotherapy-induced response (Vishnoi et al, 2020; Wang et al, 2021; Otálora-Otálora et al, 2023).