- Research article
- Open access
- Published:
Cellular reprogramming is driven by widespread rewiring of promoter-enhancer interactions
BMC Biology volume 21, Article number: 264 (2023)
Abstract
Background
Long-range interactions between promoters and cis-regulatory elements, such as enhancers, play critical roles in gene regulation. However, the role of three-dimensional (3D) chromatin structure in orchestrating changes in transcriptional regulation during direct cell reprogramming is not fully understood.
Results
Here, we performed integrated analyses of chromosomal architecture, epigenetics, and gene expression using Hi-C, promoter Capture Hi-C (PCHi-C), ChIP-seq, and RNA-seq during trans-differentiation of Pre-B cells into macrophages with a β-estradiol inducible C/EBPαER transgene. Within 1h of β-estradiol induction, C/EBPα translocated from the cytoplasm to the nucleus, binding to thousands of promoters and putative regulatory elements, resulting in the downregulation of Pre-B cell-specific genes and induction of macrophage-specific genes. Hi-C results were remarkably consistent throughout trans-differentiation, revealing only a small number of TAD boundary location changes, and A/B compartment switches despite significant changes in the expression of thousands of genes. PCHi-C revealed widespread changes in promoter-anchored loops with decreased interactions in parallel with decreased gene expression, and new and increased promoter-anchored interactions in parallel with increased expression of macrophage-specific genes.
Conclusions
Overall, our data demonstrate that C/EBPα-induced trans-differentiation involves few changes in genome architecture at the level of TADs and A/B compartments, in contrast with widespread reorganization of thousands of promoter-anchored loops in association with changes in gene expression and cell identity.
Background
Cis-regulatory elements such as enhancers and promoters of genes they control can be separated by large genomic distances [1]. However, during gene transcription, they are found in very close physical proximity [2]. Gene expression is regulated by temporal-spatial, enhancer-promoter interactions. Such long-range interactions often bypass proximal genes to exert control over specific distal target genes [1, 3]. Deciphering the principles and mechanisms underlying enhancer-promoter dynamics is essential for understanding the molecular mechanism of gene control in cell differentiation and development.
Genome-wide contact frequencies determined by Hi-C have revealed higher-order chromatin structural features. Topologically associated domains (TADs) are large megabase-size genomic regions of increased self-interaction that are partially insulated from contacts with neighboring domains [4, 5]. Various studies have reported conflicting results on TAD organization in response to environmental stimuli. Some studies show that TADs are largely invariant among different cell types and in response to environmental signals [4, 6, 7], while others have reported significant TAD re-organization in response to stimuli and between cell types [8,9,10]. Further studies have shown that transcription factors drive dynamic change of chromatin interactions between regulatory elements and genes within TADs during ESC differentiation, human fibroblast trans-differentiation, and B cell reprogramming [11,12,13].
Analysis of Hi-C heatmaps has revealed a “plaid pattern” of increased or decreased contact frequencies between domains which corresponds with their activity state. Active domains tend to contact other active domains, while inactive domains tend to contact other inactive domains, indicating that the nuclear genome is organized into two nuclear compartments referred to as A (active compartment) and B (inactive compartment) [14]. Here again, several studies have investigated compartment changes during cell differentiation in both mouse and human cell types and reported conflicting results. For example, human embryonic stem cell (ESC) differentiation to four different lineages showed substantial switching between A and B compartments [6] while during cell differentiation and development, studies show there is minimal A–B switching [15,16,17,18].
Hi-C has been instrumental in discovering and defining the principles of the above aspects of higher-order genome organization. However in the absence of billions of read-pairs, Hi-C lacks the resolution to identify significant interactions between individual promoters and their long-range regulatory elements. Promoter Capture Hi-C (PCHi-C) identifies significant interactions between individual promoters and their long-range regulatory elements in an unbiased manner [19,20,21]. PCHi-C studies have shown dynamic rewiring of enhancer-promoter interactions during human ESC differentiation to neuroectodermal cells [22], during mouse adipocyte differentiation [7], and during keratinogenesis [23] suggesting that promoter interaction changes regulate gene expression associated with cell-fate transitions. These studies are in agreement with Highly Integrative Chromatin Immunoprecipitation (HiChIP) experiments, which show chromatin reorganization during cell trans-differentiation and reprogramming [24,25,26]. However, the degree to which promoter interactions are reorganized during direct lineage conversion remains poorly understood.
Pre-B cells are precursors of B lymphocytes, which are adaptive immune cells of the lymphoid lineage [27, 28]. Macrophages are phagocytic cells of the myeloid lineage and are responsible for detecting, engulfing, and destroying pathogens [29]. The lymphoid and myeloid lineages exhibit distinct gene expression patterns maintained by lineage-restricted transcription factors [30,31,32]. The CCAAT enhancer binding protein alpha (C/EBPα) is essential for the development of the myeloid lineage [33]. An inducible Pre-B cell line (c10) with a C/EBPαER transgene has been well studied for its ability to reproducibly trans-differentiate into macrophages [34,35,36]. C/EBPα has previously been shown to bind to myeloid enhancers in Pre-B cells and regulate myeloid cell type-specific genes [33].
In this study, we used ChIP-seq, RNA-seq, Hi-C, and PCHi-C data to investigate higher-order chromatin organization and promoter interaction dynamics associated with gene expression changes during C/EBPαER-induced trans-differentiation of mouse Pre-B cells into macrophages over a 48-h period [34, 37]. Our findings show that extensive temporal rewiring of promoter-anchored interactions during the switch from a Pre-B cell-specific to a macrophage-specific transcription program constitutes the major alterations in genome organization during trans-differentiation.
Results
Transcriptional dynamics during rapid conversion of Pre-B cells into macrophages
To study the role of promoter-enhancer interactions during Pre-B cell trans-differentiation into macrophages, we employed direct lineage conversion of the C10 mouse Pre-B cell line which harbors an estradiol-inducible rat C/EBPα fused to the estrogen hormone receptor binding domain (C/EBPαER) [34, 38]. Upon addition of β-estradiol, C10 cells are converted into macrophages within 48 h. Immunofluorescence staining and western blot analysis show that C/EBPα protein is rapidly translocated from the cytoplasm to the nucleus within 1 h of β-estradiol induction and gradually decreases over 48 h (Fig. 1a, b). We verified trans-differentiation efficiency by FACS to determine the percentage of cells expressing the Pre-B cell-specific cell surface marker CD19 and macrophage-specific cell surface marker CD11b. Over 48 h of induction, the cells went from 99.5% CD19 positive to greater than 80% CD11b positive (Fig. 1c), demonstrating efficient differentiation into macrophages.
Surprisingly, we found that Histone 3 (H3) levels increased in the whole cell extracts through Pre-B trans-differentiation (Fig. 1b). To analyze gene expression dynamics during trans-differentiation, we prepared triplicate RNA-seq libraries from each of the four time points (0 h, 1 h, 12 h, and 48 h). Principal component analysis (PCA) shows three biological replicates cluster together, indicating high similarity and library quality (Additional file 1: Fig. S1a). Allele-specific analysis of H3 gene expression shows that the increased H3 protein level is due to expression of the histone variant H3.3 (H3f3a and H3f3b) genes (Fig. 1d). H3.3 is often referred to as the replacement histone as it is incorporated into transcribed chromatin throughout the cell cycle [39, 40]. Expression of canonical H3 alleles decreased dramatically after induction with CEBP/α (Fig. 1d), consistent with previous studies showing that CEBP/α arrests cell proliferation [41,42,43]. Increased levels of H3.3 may be required for the large-scale changes in gene expression that occur upon trans-differentiation. RNA-seq analysis also shows that expression of the Rat C/EBPα transgene initially increases at 1h and then decreases by 48h (Fig. 1e). The endogenous C/EBPα gene (mouse specific) was not expressed throughout the time course (Additional file 1: Fig. S1b). We identified thousands of differentially expressed genes (DEGs) comparing RNA seq data from 1h, 12h, and 48h to 0h (Additional file 1: Fig. S1c) and clustered them to reveal 5 clusters of genes with different expression profiles (Fig. 1f, Additional file 2: Table S1). Cluster 1 genes decrease in expression from 0 to 48h. Gene ontology (GO) indicates the enrichment of genes involved in B-cell-related functions: B cell activation, lymphocyte activation, and leukocyte activation (Additional file 1: Fig. S1d). A few examples of these genes are Rag1, Rag2, Hdac7, and Pax5 (Fig. 1g, Additional file 1: Fig. S1e). Cluster 2 genes increase in expression from 0h to 48h, and GO analysis shows enrichment of genes involved in the cellular response to molecule of bacterial origin, cellular response to biotic stimulus, and other macrophage-related functions (Additional file 1: Fig. S1d); some examples are Cebpb, Csf1r, Fcgr1, and Itgam (Fig. 1g and Fig.S1e). Similar to other studies, we observed downregulation of B lineage transmembrane protein gene Cd19 and B cell receptor complex signaling genes Blnk, Cd79a, Cd79b, Vpreb1, Vpreb2, and Vpreb3, while myeloid marker genes such as granulocyte collagenase 8 (Mmp8), macrophage scavenger receptor (Msr1), myeloid restricted serine protease C (Ctsc), and the myeloid cytokine-dependent chemokine 6 (Ccl6) were upregulated, (Additional file 1: Fig. S1f) [35, 44]. Together, these results show that nuclear translocation of C/EBPαER results in rapid diminution of the Pre-B cell-specific gene expression program and progressive establishment of macrophage-specific gene expression, indicating efficient conversion of Pre-B cells to macrophages.
Higher-order chromatin organization during Pre-B cell trans-differentiation
To investigate higher-order chromatin dynamics during trans-differentiation, we performed Hi-C in duplicate at the same time points. HiCUP analyses reveal > 70% of di-tags are cis-contacts greater than 10Kb, approximately 15% cis-contacts < 10 Kb, and ~ 14% trans contacts, indicating very high-quality Hi-C libraries (Additional file 1: Fig. S2a) [45]. Comparison of Hi-C matrices at 500-kb resolution showed very high similarity scores between replicates at all time points, and progressively lower similarity when comparing uninduced cells to trans-differentiation time points (Fig. 2a). Hi-C heatmaps at 100-kb and 25-kb resolution did not show gross changes in overall architecture between different time points (Additional file 1: Fig. S2b). We found approximately 3000 TADs at each time point, ranging in size from a few hundred kilobases to a few Mbs, with an average size of approximately 800kb (Additional file 1: Fig. S3a). TAD boundaries show higher reproducibility between biological replicates than between different time points (Fig. 2b), indicating a small percentage of boundaries change during trans-differentiation (Fig. 2b, c). We observed that 17,762 genes (32.08% of total genes) are located within TADs that show boundary changes (48h vs 0h). Of these genes, 937 are differentially expressed genes (DEGs) (31.36% of total DEGs). Fisher’s exact test shows that genes within TADs that show boundary changes are not enriched for DEGs (p = 0.40). GO analysis of coding genes (n = 6937) within TADs that show boundary changes shows enrichment for immune processes, such as response to chemokine, response to dsRNA, and response to pheromone (Additional file 1: Fig. S3b). These results show that a majority of TAD boundaries are conserved, consistent with previous reports of TAD conservation across different cell types [4, 6, 7]. However, we find that some boundaries do change and are associated with changes in gene expression within the flanking TADs suggesting that alteration of TAD boundaries is involved, though not a major feature driving the observed widespread changes in gene expression.
We called A/B compartments using principal component analysis (PCA). PC1 value of PCA shows a higher correlation between biological replicates than between different time points (Additional file 1: Fig. S3c). We found as expected A compartment genes have significantly higher expression levels than those in B compartments (Additional file 1: Fig. S3d). Previous reports indicate that although most compartments are stable across different cell types, some compartments do switch status during trans-differentiation in a cell-type-specific manner, reflecting cell-type-specific changes in transcription [6, 11]. We observed that a small fraction of the genome switches compartment status (either A-to-B or B-to-A) during trans-differentiation (Fig. 2d). Looking at the relationship between gene expression changes and compartment switching, we observed that of the 1688 upregulated genes from 0h to 48h, only 54 (3.2%) switched from the B compartment to the A compartment, and only 25 (2%) of 1272 downregulated genes switched from A compartment to B compartment (Fig. 2e). Greater than 90% of DEGs (either up- or downregulated) remained in stable A compartments throughout trans-differentiation (Fig. 2e and Additional file 1: Fig. S3e). Correlation between the A/B compartment switches and transcriptome clusters shows as expected A to B switches are enriched in cluster 1 which contains downregulated genes, and B to A switches are enriched in clusters 2, 3, and 4 which are composed of upregulated genes (Fig. 2f). The expression changes of genes from B to A compartments are significantly greater than those maintained in the same compartments, while the expression changes of genes from A to B compartments are notably lower (Additional file 1: Fig. S3f). GO analysis indicates that genes that switch compartments from A to B are involved in immunoglobulin production and genes that switch compartments from B to A are associated with Pre-B differentiation and macrophage-related functions, such as antifungal immune response, stimulatory C-type lectin receptor signaling pathway, and cellular response to type II interferon (Additional file 1: Fig. S3g). For example, Gbp family genes, Gbp9 and Gbp4, that are involved in response to interferon-gamma show increased expression at 48h accompanied by compartment B to A switch (Fig. 2g) [46]. Similarly, Clec4 family genes, Clec4n, Clec4d, and Clec4e, involved in macrophage defense response to other organisms also show increased expression along with the B to A compartment switch (Fig. 2g) [47]. In summary, trans-differentiation involves a limited number of changes in TAD boundaries and A/B compartments.
Promoter-anchored chromatin interactions during trans-differentiation
To investigate promoter-anchored chromatin interactions during Pre-B cell trans-differentiation, we generated Promoter Capture Hi-C (PCHi-C) libraries at each time point. PCHi-C data quality was checked using HiCUP analysis (Additional file 1: Fig S4a). Comparison of the datasets revealed high similarity scores between replicates at all time points and progressively lower similarity scores as trans-differentiation advanced (Fig. 3a). Using CHICAGO (cut-off score = 5), we identified around 130,000 significant promoter interacting regions (PIRs) at each time point (Fig. 3b). Approximately 17% are promoter to promoter (bait to bait), and 83.4% are promoter to other genomic fragments, (Additional file 1: Fig. S4b). Promoter-anchored interactions across all time points show similar distributions over linear genomic distance (Fig. 3b) with median distances of 163 ~ 184kb.
Cis-regulatory elements, such as enhancers, far outnumber genes in the mammalian genome [48, 49]. In many cases, a gene will be contacted by multiple enhancers, which can be located at great genomic distances [1, 50]. We examined the relationship between the number of significant interactions per gene promoter and gene expression level. We found that genes possessing a larger number of interactions tend to have significantly higher gene expression (Fig. 3c) in agreement with previous studies [50,51,52].
Functionally active cis-regulatory elements are often enriched in specific histone modifications [53]. H3K27ac, H3K4me1, and p300 are often associated with active enhancers, whereas H3K27me3 modifications are generally associated with Polycomb repressed chromatin, bivalent domains, and poised enhancers [54,55,56]. To characterize PIRs during Pre-B cell trans-differentiation, we analyzed publicly available ChIP-seq data for H3K27ac, H3K4me1, H3K27me3, and p300 peaks [37] at PIRs compared to distance-matched, random control regions at all four time points. We found that PIRs are significantly enriched for active enhancer marks H3K27Ac, H3K4Me1, and p300 (Fig. 3d), consistent with active transcription of the contacted promoters. We also found that PIRs are significantly enriched with H3K27me3 peaks (Fig. 3d). To investigate the role of PIRs enriched in H3K27ac and H3K27me3 in gene expression, we classified promoters into four groups: Those with at least one PIR marked by H3K27ac, at least one PIR marked by H3K27me3, both, or no marks. Promoters contacted by PIRs with H3K27ac show significantly higher gene expression than the other three groups. Promoters contacted by PIRs with H3K27me3 show significantly lower gene expression than the other three groups (Additional file 1: Fig. S4c). We classified promoters into four groups based on the number of their PIRs marked by H3K27ac and compared their expression level. The results show that promoters connected to more than one PIR marked by H3K27ac peaks have significantly higher gene expression indicating a dose-dependent effect on transcription (Fig. 3e). To quantify the effect of H3K27me3, we focused on promoters contacting PIRs marked by H3K27me3 and not marked by H3K27ac. These promoters were then categorized into four groups based on the number of PIRs marked by H3K27me3. The results show that promoters that connected more than one PIR marked by H3K27me3 have a modest but significantly lower gene expression level than promoters that only have one PIR marked by H3K27me3, suggesting a dose-dependent effect of repressor-promoter interactions on gene expression (Additional file 1: Fig. S4d). These results demonstrate that active and repressive histone modifications at PIRs are associated with active and repressed transcription (respectively) of the promoters they contact, as expected.
Dynamic rewiring of interactions during trans-differentiation
We measured significant differences in promoter interactions during trans-differentiation using ChicDiff [57] and found that the number of differential interactions increased dramatically during trans-differentiation. Compared to uninduced cells, we found 47 differential interactions involving 18 promoters at 1h, 862 differential interactions at 225 promoters at 12h, and 5847 differential interactions at 1432 promoters at 48h (Fig. 4a), Additional file 3: Table S2). To investigate the potential effect of differential interactions on gene expression, we classified promoters into two groups based on contact frequency (read-pair depth): increased (gain of contacts at 1h, 12h, or 48h compared to 0h) and decreased (decreased contact frequency between 1h, 12h, and 48h compared to 0h). These data show that contact frequency is positively correlated with gene expression. Promoters that gained contacts show significantly upregulated gene expression while promoters that lost contacts show significantly downregulated gene expression (Fig. 4b). To investigate the role of histone modifications on gene expression by chromatin looping, we classified the PIRs of significant differential promoter interactions into three categories based on overlap with different histone marks and quantified their target gene expression change. We found gene promoters that gained chromatin interactions with enhancer marks (either H3K27ac or H3K4me1) show significantly upregulated gene expression compared to gene promoters that only gain contact with repressor marks (H3K27me3) (Fig. 4c). Gene promoters that gained contact with repressor marks (H3K27me3) such as Maf1 and Sspb1 show downregulated gene expression (Additional file 1: Fig. S5a, Additional file 1: Fig. S5b). We also investigated the additive effect of gaining new enhancer contacts on gene expression changes by quantifying the number of gained PIRs with H3K27ac or H3K4me1 peaks and evaluated the corresponding gene expression changes at 48h compared to 0h. Promoters gaining interactions with multiple PIRs with enhancer marks exhibited significantly higher gene expression changes compared to those with interactions involving only one enhancer mark (Fig. 4d), Additional file 1: Fig. S5c). Taken together, these results suggest that gene expression changes during trans-differentiation are driven in part through widespread changes in promoter interactions with functional regulatory elements.
We clustered differential interactions to further investigate the dynamic rewiring of promoter interactions (Fig. 4e), Additional file 3: Table S2). The most dramatic change in differential interactions occurs between 12h and 48h. Gene ontology (GO) enrichment analysis shows that these genes are associated with macrophage-related functions, including response to chemokine, interferon-gamma, and cytokine, as well as cell migration and motility (Fig. 4f). For example, Ccl7 and Ccl12 gain interactions at 48h and have macrophage-related functions (Fig. 4g) [58, 59]. Interestingly, GO shows genes that rapidly lose interactions during trans-differentiation are involved in regulation of chromosome organization, cell cycle, and negative regulation of transcription, suggesting alterations in chromosome organization and cell cycle changes occur along with changes in gene expression (Fig. 4f). Suz12 is an example of a gene that progressively loses interactions after induction (Fig. 4g). We identified gene promoters that lost interactions at 48h, and GO analysis shows enrichment of genes involved in type I interferon production and B cell activation (Additional file 1: Fig. S6a). Some examples of Pre-B cell-specific genes are Pax5 and Clcf1 (Additional file 1: Fig. S6b). These results indicate that gained and lost promoter-anchored interactions are a major contributor to cell fate transition.
C/EBPα binds to interacting regions and modulates gene expression during Pre-B cell trans-differentiation
The C/EBPα transcription factor is a master regulator of macrophage identity [33]. We investigated how C/EBPα binding to interacting regions affects transcriptional changes. We identified 23,434 C/EBPα binding peaks 1h after induction, and a gradual decrease at 12h and 48h (Additional file 1: Fig. S7a). At 0h, C/EBPα binds only to ~ 4% of significant interactions (either promoters or PIRs). After induction, C/EBPα binds to ~ 29% of significant interactions at 1h, ~ 18% at 12h, and ~ 13% at 48h (Fig. 5a). We observed that PIRs have a significantly higher level of C/EBPα occupancy compared to random DNA controls (Fig. 5b). To investigate the association between differential C/EBPα binding and gene expression, we identified time-point unique and common C/EBPα binding sites at 1h, 12h, and 48h, compared to 0h (Additional file 1: Fig. S5b). We analyzed the expression fold change of genes with differential C/EBPα occupancy at their promoters. At 1h, 12h, and 48h post-C/EBPα induction, C/EBPα-occupied DEG promoters showed significantly increased expression (Fig. 5c). At 12h and 48h, increased C/EBPα occupancy at PIRs also correspond to increased expression levels of the DEGs they contact (Fig. 5d). For example, Btg1 located on chr10, involved in cell growth and differentiation [60, 61], gradually gains contacts both upstream and downstream as trans-differentiation progresses (Fig. 5e). Several gained PIRs over a 490-kb region upstream of the Btg1 promoter show increased C/EBPα binding and gain of H3K27ac (Fig. 5e). The Btg1 promoter itself is not occupied by C/EBPα, suggesting that the several C/EBPα-dependent enhancers in the region modulate Btg1 expression through interaction with the Btg1 promoter. Taken together, these results show numerous direct targets of C/EBPα binding to promoters and PIRs of differentially expressed genes, suggesting that changes in promoter interaction loops are a major driver of gene expression changes during Pre-B trans differentiation into macrophages.
Discussion
Consistent with previous studies [34, 35], we observed efficient trans-differentiation of the C/EBPαER-expressing B cell line C10, providing us with a robust and reproducible system to investigate transcription factors driven chromosome reorganization during direct reprogramming from Pre-B cells to macrophages. Nuclear translocation of C/EBPα in Pre-B cells results in differential expression of thousands of genes. Transcription of B-cell-related genes is switched off and genes associated with macrophage functions and identity are switched on over the 48h period tested. We show that C/EBPα binding to thousands of sites including promoters and PIRs results in a rapid, significant change in the promoter interactome, with C/EBPα binding to 30% of significant promoter interactions genome-wide within 1h of induction. Genes with induced C/EBPα binding to their promoters are the first to respond showing significantly increased expression within 1h of induction. Genes whose promoters gain interactions with distal C/EBPα binding sites show delayed induction showing significantly increased expression at 12h. Our data also suggest widespread indirect effects of C/EBPα at other DEGs are associated with changes in enhancer-promoter interactions.
Altered gene expression during C/EBPα-induced trans-differentiation could not be explained by widespread changes in TAD or A/B compartment organization. We observed very high conservation of TAD boundary locations genome-wide with only a small percentage that shift location in association with the cell fate transition. Similarly, although some genes associated with macrophage functions, such as Gbp9 and Gbp4, switch from the B compartment to the A compartment during trans-differentiation, the vast majority of DEGs (94%) did not change compartment status over the 48h period. A recent study that used C/EBPα to trans-differentiate a human B cell line to macrophage also shows conservation of TAD boundaries, and observed only a small percentage of dynamic compartments associated with gene expression changes [62]. It is important to point out that these cells retain macrophage morphology, functions, and gene expression even after removal of estradiol indicating a stable change of cell fate [34, 63].
Our results suggest that the rewiring of promoter-anchored interactions plays a major role in controlling differential gene expression and cell fate transition compared to the relatively minor role played by changes in higher-order organization of TADs and A/B compartments. We show that PIRs are significantly enriched in histone modifications associated with gene regulatory activity. Genes whose promoters are contacted by PIRs with the active enhancer mark, H3K27ac, are highly expressed, while promoters contacted by the repressive H3K27me3 modification are not. Thus, transcription factor-mediated differential loop formation between promoters and their distal regulatory elements is the major driver of the cell-type-specific, promoter interactomes [50] and is the likely basis for tissue-specific gene expression profiles.
Conclusions
Our integrated analyses of gene expression, histone modifications, transcription factor binding, and various levels of three-dimensional chromosome conformation during trans-differentiation have revealed that widespread changes in gene expression that characterize the transition from a Pre-B cell to a macrophage phenotype involve surprisingly few changes in TAD structure or A/B compartment status. Instead, we found reorganization of thousands of promoter-anchored loops, in conjunction with transcription factor binding and histone modification dynamics, occurring at thousands of differentially expressed genes suggesting that changes in chromatin interactions at the sub-TAD level are a major mechanism of differential regulation of gene expression that drives cell state transitions.
Methods
Cell culture, reprogramming induction, and flow cytometry
C10 cells were cultured in RPMI 1640 (ThermoFisher, 11875–093) medium supplemented with 10% heat-inactivated FBS (ThermoFisher, A3840202), 1% Penicillin–Streptomycin (ThermoFisher, 15140122), and 50 μM 2-Mercaptoethanol (ThermoFisher, 31350010). For C10 cells trans-differentiating into macrophages, 5.0 × 106 C10 cells were plated in 10-cm cell culture dishes, followed by adding100 nM of β-estradiol (Sigma, E2758-250MG), 10 nM IL-3 (PEPRO TECH, 213–13), and CSF-1 (PEPRO TECH, 315–02). The cell culture medium was replenished every 24h. The Pre-B cell trans-differentiation efficiency was measured by flow cytometry analysis and qPCR gene expression analysis of Pre-B cell-specific genes and macrophage-specific genes. For flow cytometry analysis, C10 cells at various differentiation stages were stained using PE anti-mouse/human CD11b (BioLegend, 101208) antibody and APC anti-mouse CD19 antibody (BioLegend, 152410). The percentage of Pre-B cells and macrophage population was analyzed on FACSCanto II analyzer (BD Biosciences). Flow cytometry data was analyzed by FlowJo (TreeStar, Inc) software. PE Rat IgG2b (BioLegend, 400508) and APC Rat IgG2a (BioLegend, 400511) antibodies were used as isotype controls to exclude the background fluorescence.
Western blot and immunofluorescence staining
For western blots, cells were lysed by CelLytic buffer (Millipore Sigma, C3228). The protein concentration was quantified by the Micro-BCA assay reagent (ThermoFisher Scientific, 23235). A 10 μg protein from each sample was loaded into 4–12% SDS–polyacrylamide (ThermoFisher Scientific, NP0321Box) for electrophoresis. The following antibodies were used for Western blot experiments: Anti-C/EBPα (Santa Cruz Biotechnology, AF2018), Anti-H3 (Cell Signaling Technology, 9715S), Anti-β-Tubulin (Cell Signaling Technology, 2128S). For immunofluorescence, C10 cells were first fixed in 4% PFA for 5 min, washed using cold PBS twice, followed by attaching cells to fibronectin and pre-coating the glass slides. Next, the cells were treated using blocking buffer (5% normal goat serum and 0.3% Triton X-100 in 1X PBS) for an hour at room temperature, followed by primary antibody treatment overnight at 4 °C and secondary antibody treatment at room temperature for an hour. The primary and secondary antibodies were diluted in antibody dilution buffer (1% BSA and 0.3% Triton X-100 in 1X PBS), and DAPI (Sigma Millipore, D9542) was used for nuclear staining. Anti-C/EBPα (Santa Cruz Biotechnology, AF2018) and Goat anti-Mouse IgG (H + L) Alexa Fluor 488 (ThermoFisher, A-1100A) antibodies were used for immunofluorescence staining. All the immunofluorescence images were taken by LSM 710 confocal microscope (Zeiss) and processed using Zen black (Zeiss).
RNA extraction, qPCR assay of gene expression, RNA-seq libraries preparation, and sequencing
RNA was prepared from C10 cells at various trans-differentiation stages using RNeasy mini kit (Qiagen, Cat#74104) followed by on-column DNase I treatment. Total RNA (1 μg) was used for reverse transcription to generate cDNA by using iScript™ Reverse Transcription Supermix for RT-qPCR (BioRad, Cat#1708840), followed by the real-time PCR amplification using Power SYBR Green PCR Master Mix (Thermofisher SCIENTIFIC, Catalog number 4367659). The real-time PCR was performed with QuantStudio 7 Flex Real-Time PCR System (Thermofisher SCIENTIFIC) for 40 cycles, 95°C for 15 sec, and 60°C for 1 min. Gene expression was normalized to the housekeeping gene Gapdh using the 2−ΔΔCT method.
For RNA-seq libraries preparation, 1μg of total RNA (RIN > 7) was used for rRNA depletion by using NEB rRNA Depletion Kit (E6310). Next, the NEBNext® Ultra™ II RNA Library Prep Kit for Illumina (E7770S) and NEBNext NEBNext® Multiplex Oligos for Illumina® (E7335) were used for library construction. RNA-seq library concentration was quantified by using the TapeStation and KAPA Library QANT Kit (Roche, 07960336001), followed by single-end (1 × 100bp) sequencing on the Novaseq S2 flow cell. The sequencing depth is described in Additional file 2: Table S1.
Hi-C library and Promoter Capture Hi-C (PCHi-C) library preparation, sequencing
Hi-C and Promoter Capture Hi-C (PCHi-C) library preparation was described in our previous publication with minor modifications [20]. Ten million cells were used for each Hi-C library preparation. Briefly, we fixed cells using 2% formaldehyde for 10 min, followed by quenching in 0.125 M glycine for 5 min on ice. The fixed cells were washed twice in cold PBS and centrifuged at 760 g for 5 min at 4°C. Next, the cells were lysed in 20 mL cold lysis buffer (10 mM Tris–HCl pH 8, 0.2% IGEPAL CA-630, 10 mM NaCl, and one tablet protease inhibitor cocktail) on ice for 30 min, followed by centrifuging at 760 g for 5 min to remove the supernatant. The pellet was washed in 1.25X NEB buffer2 (NEB, B7002S) and resuspended in 358 μL 1.25X NEB buffer2, followed by adding 11 μL of 10% SDS and incubating at 37 °C with shaking at 950 rpm for 30 min. Lastly, the SDS was quenched by adding 75μL of 10% Triton X-100 and incubating at 37 °C for 15 min with shaking. To digest chromatin using HindIII restriction enzyme, we added 12 μL of 100 U/μL HindIII (NEB, R0104M) to each reaction, followed by shaking at 37 °C and 950 rpm overnight, followed by adding 5 μL of 100 U/μL HindIII (500 units in total) per reaction at 37 °C for 2 h on the following morning. To repair the digested overhangs, we added 6.1 μL 10X NEB buffer2, 25 μL molecular-grade water, 15.3 μL of 1 mM biotin-14-dATP (Jena Bioscience, NU-835-Bio14-L), 1.56 μL of 10 mM dCTP, 1.56 μL of 10 mM dGTP, 1.56 μL of 10 mM dTTP, and 10.2 μL of 5U/μL DNA polymerase I, Large (Klenow) Fragment (NEB, M0210L) and incubated the reaction at 37 °C for 1 h. In-nucleus ligation mix was prepared by adding 102μL 10X T4 DNA ligase buffer, 10.2 μL BSA (NEB, B9000S), 350.9 μL molecular-grade water, and 25.5 μL 1U/ μL T4 DNA ligase (ThermoFisher Scientific, 15224017), followed by incubation at 16 °C for 4 h in a thermomixer (shaking at 700 rpm for 5 s in every 2 min). Next, the nuclei were pelleted at 2500 g for 5 min and resuspended in 300μL cross-link reversal buffer (10 mM Tris–HCl, 0.5 M NaCl, 1% SDS), followed by adding 5 μL 10 mg/ml RNase A (ThermoFisher Scientific, EN0531) at 37 °C for 30 min and 20 μL of 20 mg/ml Proteinase K (Gold Bio, P-480-SL2) at 55 °C for 1 h and 68 °C overnight. The genomic DNA was extracted using cold ethanol and sodium acetate (pH 5.2) and resuspended in 130μL of 10 mM Tris–HCl. DNA was transferred to a microTUBE AFA fiber Pre-Slit Snap-Cap (Covaris, 520045) and fragmented using a ME220 Covaris sonicator (Peak Incident Power 50W, Duty Factor 20%, Cycles per Burst 200, Treatment time 80 s). The sonicated DNA was transferred to a 1.5-ml tube, and the total volume was brought to 200 μL by adding 70 μL of molecular-grade water. DNA size was double size selected using AMPure XP beads (BECKMAN COULTER, A63881). First, 120 μL of beads were added to each reaction (ratio of AMPure beads to DNA: 0.6 to 1) and incubated for 5 min. The clear supernatant was transferred to a fresh tube through a magnetic separation stand. Next, 30 μl of fresh AMPure XP beads were added to the clear supernatant and incubated for 5 min at room temperature. Beads were separated on a magnetic separation stand and washed twice using 800 μL of 70% ethanol and dried at 37 °C, followed by eluting DNA in 300 μL 1X Tris buffer. DNA was used for library preparation with a normal size distribution between 300 and 500 bp and was visualized by TapeStation.
Procedures for Biotin/Streptavidin pull-down of the Hi-C ligation products, end repair, removal of biotin at the non-ligated DNA ends, adaptors ligation, Hi-C libraries amplification, and Hi-C size selection were described in our recent publication [64]. The Hi-C library size distribution was evaluated by using TapeStation, and the Hi-C library concentration was estimated by using the Qubit 4 Fluorometer and KAPA Library QANT Kit (Roche, E7770S). Procedures and baited capture system design for Promoter Capture Hi-C libraries preparation were described in our previous publication [19, 20]. 39,021 biotinylated RNA bait target 22,225 annotated gene promoters that include protein-coding, non-coding(lincRNA), antisense, snRNA, miRNA, or snoRNA. Two unique 120 bp capture probes were designed close to the ends of each restriction fragment (one to each end) containing a transcription start site. In rare instances, a unique sequence could not be found for one end, or even more rare, for both ends. For each reaction, 500 ng of dried Hi-C DNA was used for library preparation. SureSelectXT Custom 3–5.9 Mb (Agilent, 5190–4831) and SSEL TE Reagent Kit, ILM PE FULL Adaptor (Agilent, 931108) were used for enriching the Hi-C fragment that is associated with promoters. The procedures for Promoter Capture Hi-C libraries amplification were described in detail in our recent publication [64]. The Promoter Capture Hi-C library size distribution was evaluated by using TapeStation, and the library concentration was quantified by using the Qubit 4 Fluorometer and KAPA Library QANT Kit (Roche, E7770S). Hi-C and Promoter Capture Hi-C libraries were paired-end sequenced (2 × 50 bp) on Novaseq S2 flow cells. The sequencing depth is described in Additional file 1: Table S1.
RNA-seq data processing
RNA-seq libraries were prepared for each time point. We included three biological and four technical replicates for each time point. Technical replicates were merged, and Trim Galore (0.6.2) pipeline was used to remove adapters and short reads from the data (https://github.com/FelixKrueger/TrimGalore). FastQC program was used to check the quality of the data (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Trimmed reads were mapped to mouse reference genome (mm10) using Hisat2 (2.1.0) with default parameters [65]. Samtools (1.10) were used to convert the sam files to bam files [66]. GTF file (GRCm38.96) was used to annotate the mapped reads using featureCounts tool [67]. Exon was chosen as the feature gene type with a default minimum overlap of 1 bp. The genes with the sum of read number in all samples less than 2 were removed from further analysis. Differentially expressed genes (DEGs) were identified using R package DEseq2 by keeping filtering thresholds of Padj < 0.05 and fold change > 2 [68]. Sample distances were assessed and visualized by the principal components analysis (PCA) function in R. TPM (transcripts per million) values of each gene were obtained using RSEM [69]. Fuzzy c-means clustering was performed on DEGs using R [70]. GO analysis was performed on WebGestalt toolkit with “over-representation analysis” as the method of interest and “biological process” as the functional database [71]. SNPsplit was used to analyze the expression of the rat Cebpa transgene and endogenous mouse Cebpa gene [72].
Hi-C data analysis
Hi-C data were processed using HiCUP (v0.7.2) pipeline [73] with mm10 as the reference genome. Hi-C interaction matrices’ reproducibility of the biological replicates was assessed by Hi-CRep (v1.8.0) using the stratum-adjusted correlation coefficient method [74]. Hi-CUP bam files were converted to Hi-C files using Hi-Cup2homer and tagDir2Hi-CFile.pl script in HiCUP and Homer package(http://homer.ucsd.edu/homer/). JuiceBox was used to visualize Hi-C matrix heatmap and TAD structure [75]. Homer package was used to call the A/B compartment at 50-kb resolution [76]. The correlation coefficient of PC1 value was calculated to assess the reproducibility of biological replicates. A/B compartment switches were estimated based on PC1 values. To calculate the correlation between DEGs and A/B compartment switch, A/B switch was classified into four categories: A to B, B to A, A stable, and B stable. DEGs were assigned to each category based on overlap between the start coordinates of DEGs and compartment bin. The A/B compartment was visualized using Integrative Genomics Viewer (IGV) [77]. Hi-C files were converted to h5 format at 50-kb resolution, and KR correction was performed using Hi-CExplorer (v3.3) tool [78]. TADs and TAD boundaries were defined by Hi-CExplorer at 50-kb resolution with parameters (Hi-CFindTADs –correctForMultipleTesting fdr –thresholdComparisons 0.01 –delta 0.01) [78].
PCHi-C interaction analysis
PCHi-C mouse raw reads were processed by the HiCUP pipeline and were mapped to mm10 genome [73]. Hi-CRep was applied to assess the PCHi-C interaction matrices’ reproducibility of the biological replicates [74]. CHICAGO (v1.12.0) pipeline was used to call the significant interactions with CHiCAGO score cut-off 5 [79]. PCHi-C contacts were visualized using WashU Browser (https://epigenomegateway.wustl.edu/) [80]. To estimate the correlation between gene expression and the PCHi-C interactions, the expression level was classified into 4 categories: TPM = 0, TPM = 0–1, TPM = 1–5, and TPM > 5. The enrichment of markers (H3K27ac, H3K4me1, p300, and H3K27me3) on promoter interacting regions (PIRs) was tested by shuffling the distance-matched fragments around the genome as control (peakEnrichemnt4Features function in CHICAGO) [79].
BEDTools (v2.26.0) “intersect” function was used to find PIR regions that overlap with H3K27ac peaks or H3K27me3 [81]. To investigate the association between promoter expression level and its PIRs, promoters were categorized into four groups based on whether at least one PIR overlapped with H3K27ac or H3K27me3 peaks and compared the TPM values of the respective promoters. The contacts showing significant changes at different time points were identified by running Chicdiff (v0.5) package with default setting [57]. When assessing the association between alteration of interactions and regulation of gene expression, promoters were divided into two groups based on whether their contact frequency significantly increased or significantly decreased. Pheatmap function in R was used to cluster differential interactions with a score greater than or equal to 5 in at least one time point based on the Euclidean distance for CHICAGO scores across all time points [70].
ChIP-seq data analysis
ChIP-seq data C/EBPα along with histone modification ChIP data were downloaded from GEO [37]. Raw reads were mapped end-to-end to the mm10 reference genome using Bowtie2 [82]. SAM format files were converted to BAM format and were sorted by Samtools [66]. The findPeaks function of HOMER package was used to call histone modification peaks by setting peak size 1000 bp and the remaining parameters were set to the default values [76]. Peak calling for C/EBPα data was performed using MACS2 (v2.1.1) (minimum FDR cutoff 0.01) [83]. To visualize the peak, peak files were loaded to IGV [77]. The enrichment analysis of C/EBPα binding peaks on PIRs was performed by running CHICAGO peakEnrichmen4Features [79]. Further, to quantitatively compare differential binding sites of C/EBPα at different time points after induction, MAnorm (v1.3.0) package was used with default settings [84]. p-value < 0.05 and fold change > 1.5 were used as threshold to identify induced C/EBPα binding sites. TSS of all genes were downloaded from Ensembl biomart (v100) [85]. Bedtools intersect was used to find at least 1-bp overlap between differential C/EBPα binding sites and TSSs (a total of 2–1 kb upstream and 1 kb downstream of TSS) or PIRs [81].
Statistical analysis
Unpaired t-test, Mann–Whitney U statistical test, and chi-squared test were performed using GraphPad Prism 8 and R. p-value < 0.05 was considered as significant unless specified and indicated as * p < 0.05, ** p < 0.01 and *** p < 0.001.
Availability of data and materials
Histone modifications ChIP-seq data is available at GEO with accession number GSE53173 [86]. C/EBPα ChIP-seq data is available at GEO with accession number GSE53362 [87]. RNA-seq, Hi-C, and Promoter Capture Hi-C data of Pre-B cell differentiation were deposited in Gene Expression Omnibus: GSE221737 [88].
Abbreviations
- PCHi-C:
-
Promoter capture Hi-C
- TADs:
-
Topological associated domains
- ESC:
-
Embryonic stem cell
- HiChIP:
-
Highly integrative chromatin immunoprecipitation
- H3:
-
Histone 3
- PCA:
-
Principal component analysis
- DEGs:
-
Differentially expressed genes
- GO:
-
Gene ontology
References
Schoenfelder S, Fraser P. Long-range enhancer–promoter contacts in gene expression control. Nat Rev Genet. 2019;20:437–55.
Carter D, Chakalova L, Osborne CS, Dai Y, Fraser P. Long-range chromatin regulatory interactions in vivo. Nat Genet. 2002;32:623–6.
Pennacchio LA, Bickmore W, Dean A, Nobrega MA, Bejerano G. Enhancers: five essential questions. Nat Rev Genet. 2013;14:288–95.
Dixon JR, Selvaraj S, Yue F, Kim A, Li Y, Shen Y, et al. Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature. 2012;485:376–80.
Nora EP, Lajoie BR, Schulz EG, Giorgetti L, Okamoto I, Servant N, et al. Spatial partitioning of the regulatory landscape of the X-inactivation centre. Nature. 2012;485:381–5.
Dixon JR, Jung I, Selvaraj S, Shen Y, Antosiewicz-Bourget JE, Lee AY, et al. Chromatin architecture reorganization during stem cell differentiation. Nature. 2015;518:331–6.
Siersbæk R, Madsen JGS, Javierre BM, Nielsen R, Bagge EK, Cairns J, et al. Dynamic rewiring of promoter-anchored chromatin loops during adipocyte differentiation. Mol Cell. 2017;66:420-435.e5.
Krijger PHL, Di Stefano B, de Wit E, Limone F, van Oevelen C, de Laat W, et al. Cell-of-origin-specific 3D genome structure acquired during somatic cell reprogramming. Cell Stem Cell. 2016;18:597–610.
Winick-Ng W, Kukalev A, Harabula I, Zea-Redondo L, Szabó D, Meijer M, et al. Cell-type specialization is encoded by specific chromatin topologies. Nature. 2021;599:684–91.
Li L, Lyu X, Hou C, Takenaka N, Nguyen HQ, Ong CT, et al. Widespread rearrangement of 3D chromatin organization underlies polycomb-mediated stress-induced silencing. Mol Cell. 2015;58:216–31.
Stadhouders R, Vidal E, Serra F, Di Stefano B, Le Dily F, Quilez J, et al. Transcription factors orchestrate dynamic interplay between genome topology and gene regulation during cell reprogramming. Nat Genet. 2018;50:238–49.
Dall’Agnese A, Caputo L, Nicoletti C, di Iulio J, Schmitt A, Gatto S, et al. Transcription factor-directed re-wiring of chromatin architecture for somatic cell nuclear reprogramming toward trans-differentiation. Mol Cell. 2019;76:453-472.e8.
Sun F, Chronis C, Kronenberg M, Chen X-F, Su T, Lay FD, et al. Promoter-enhancer communication occurs primarily within insulated neighborhoods. Mol Cell. 2019;73:250-263.e5.
Lieberman-Aiden E, van Berkum NL, Williams L, Imakaev M, Ragoczy T, Telling A, et al. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science. 2009;326:289–93.
Chan WF, Coughlan HD, Zhou JHS, Keenan CR, Bediaga NG, Hodgkin PD, et al. Pre-mitotic genome re-organisation bookends the B cell differentiation process. Nat Commun. 2021;12:1344.
Johanson TM, Lun ATL, Coughlan HD, Tan T, Smyth GK, Nutt SL, et al. Transcription-factor-mediated supervision of global genome architecture maintains B cell identity. Nat Immunol. 2018;19:1257–64.
Won H, de la Torre-Ubieta L, Stein JL, Parikshak NN, Huang J, Opland CK, et al. Chromosome conformation elucidates regulatory relationships in developing human brain. Nature. 2016;538:523–7.
Kurotaki D, Kikuchi K, Cui K, Kawase W, Saeki K, Fukumoto J, et al. Chromatin structure undergoes global and local reorganization during murine dendritic cell development and activation. Proc Natl Acad Sci. 2022;119:e2207009119.
Schoenfelder S, Furlan-Magaril M, Mifsud B, Tavares-Cadete F, Sugar R, Javierre BM, et al. The pluripotent regulatory circuitry connecting promoters to their long-range interacting elements. Genome Res. 2015;25:582–97.
Schoenfelder S, Javierre B-M, Furlan-Magaril M, Wingett SW, Fraser P. Promoter capture Hi-C: high-resolution, genome-wide profiling of promoter interactions. J Vis Exp. 2018. https://doi.org/10.3791/57320.
Sahlén P, Abdullayev I, Ramsköld D, Matskova L, Rilakovic N, Lötstedt B, et al. Genome-wide mapping of promoter-anchored interactions with close to single-enhancer resolution. Genome Biol. 2015;16:156.
Freire-Pritchett P, Schoenfelder S, Várnai C, Wingett SW, Cairns J, Collier AJ, et al. Global reorganisation of cis-regulatory units upon lineage commitment of human embryonic stem cells. Elife. 2017;6:e21926.
Rubin AJ, Barajas BC, Furlan-Magaril M, Lopez-Pajares V, Mumbach MR, Howard I, et al. Lineage-specific dynamic and pre-established enhancer-promoter contacts cooperate in terminal differentiation. Nat Genet. 2017;49:1522–8.
Di Giammartino DC, Kloetgen A, Polyzos A, Liu Y, Kim D, Murphy D, et al. KLF4 is involved in the organization and regulation of pluripotency-associated three-dimensional enhancer networks. Nat Cell Biol. 2019;21:1179–90.
Chen L, Cao W, Aita R, Aldea D, Flores J, Gao N, et al. Three-dimensional interactions between enhancers and promoters during intestinal differentiation depend upon HNF4. Cell Rep. 2021;34:108679.
Ren R, Fan Y, Peng Z, Wang S, Jiang Y, Fu L, et al. Characterization and perturbation of CTCF-mediated chromatin interactions for enhancing myogenic transdifferentiation. Cell Rep. 2022;40: 111206.
Gathings WE, Lawton AR, Cooper MD. Immunofluorescent studies of the development of pre-B cells, B lymphocytes and immunoglobulin isotype diversity in humans. Eur J Immunol. 1977;7:804–10.
Raff MC, Megson M, Owen JJ, Cooper MD. Early production of intracellular IgM by B-lymphocyte precursors in mouse. Nature. 1976;259:224–6.
Wynn TA, Chawla A, Pollard JW. Macrophage biology in development, homeostasis and disease. Nature. 2013;496:445–55.
Scott EW, Simon MC, Anastasi J, Singh H. Requirement of transcription factor PU.1 in the development of multiple hematopoietic lineages. Science. 1994;265:1573–7.
Lukin K, Fields S, Hartley J, Hagman J. Early B cell factor: regulator of B lineage specification and commitment. Semin Immunol. 2008;20:221–7.
Kee BL, Quong MW, Murre C. E2A proteins: essential regulators at multiple stages of B-cell development. Immunol Rev. 2000;175:138–49.
Suh HC, Gooya J, Renn K, Friedman AD, Johnson PF, Keller JR. C/EBPalpha determines hematopoietic cell fate in multipotential progenitor cells by inhibiting erythroid differentiation and inducing myeloid differentiation. Blood. 2006;107:4308–16.
Bussmann LH, Schubert A, Vu Manh TP, De Andres L, Desbordes SC, Parra M, et al. A robust and highly efficient immune cell reprogramming system. Cell Stem Cell. 2009;5:554–66.
Xie H, Ye M, Feng R, Graf T. Stepwise reprogramming of B cells into macrophages. Cell. 2004;117:663–76.
Di Tullio A, Vu Manh TP, Schubert A, Castellano G, Månsson R, Graf T. CCAAT/enhancer binding protein alpha (C/EBP(alpha))-induced transdifferentiation of pre-B cells into macrophages involves no overt retrodifferentiation. Proc Natl Acad Sci U S A. 2011;108:17016–21.
van Oevelen C, Collombet S, Vicent G, Hoogenkamp M, Lepoivre C, Badeaux A, et al. C/EBPα activates pre-existing and de novo macrophage enhancers during induced pre-B cell transdifferentiation and myelopoiesis. Stem Cell Reports. 2015;5:232–47.
Landschulz WH, Johnson PF, Adashi EY, Graves BJ, McKnight SL. Isolation of a recombinant copy of the gene encoding C/EBP. Genes Dev. 1988;2:786–800.
Shi L, Wen H, Shi X. The histone variant H3.3 in transcriptional regulation and human disease. J Mol Biol. 2017;429:1934–45.
Krimer DB, Cheng G, Skoultchi AI. Induction of H3.3 replacement histone mRNAs during the precommitment period of murine erythroleukemia cell differentiation. Nucleic Acids Res. 1993;21:2873–9.
Wang H, Iakova P, Wilde M, Welm A, Goode T, Roesler WJ, et al. C/EBPalpha arrests cell proliferation through direct inhibition of Cdk2 and Cdk4. Mol Cell. 2001;8:817–28.
Timchenko NA, Harris TE, Wilde M, Bilyeu TA, Burgess-Beusse BL, Finegold MJ, et al. CCAAT/enhancer binding protein α regulates p21 protein and hepatocyte proliferation in newborn mice. Mol Cell Biol. 1997;17:7353–61.
Slomiany BA, D’Arigo KL, Kelly MM, Kurtz DT. C/EBPα inhibits cell growth via direct repression of E2F-DP-mediated transcription. Mol Cell Biol. 2000;20:5986–97.
Francesconi M, Di Stefano B, Berenguer C, de Andrés-Aguayo L, Plana-Carmona M, Mendez-Lago M, et al. Single cell RNA-seq identifies the origins of heterogeneity in efficient cell transdifferentiation and reprogramming. Elife. 2019;8:e41627.
Nagano T, Várnai C, Schoenfelder S, Javierre B-M, Wingett SW, Fraser P. Comparison of Hi-C results using in-solution versus in-nucleus ligation. Genome Biol. 2015;16:175.
Haque M, Siegel RJ, Fox DA, Ahmed S. Interferon-stimulated GTPases in autoimmune and inflammatory diseases: promising role for the guanylate-binding protein (GBP) family. Rheumatology (Oxford). 2021;60:494–506.
Miyake Y, Toyonaga K, Mori D, Kakuta S, Hoshino Y, Oyamada A, et al. C-type lectin MCL is an FcRγ-coupled receptor that mediates the adjuvanticity of mycobacterial cord factor. Immunity. 2013;38:1050–62.
ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;489:57–74.
Shen Y, Yue F, McCleary DF, Ye Z, Edsall L, Kuan S, et al. A map of the cis-regulatory sequences in the mouse genome. Nature. 2012;488:116–20.
Javierre BM, Sewitz S, Cairns J, Wingett SW, Várnai C, Thiecke MJ, et al. Lineage-specific genome architecture links enhancers and non-coding disease variants to target gene promoters. Cell. 2016;167:1369-1384.e19.
Song M, Yang X, Ren X, Maliskova L, Li B, Jones IR, et al. Mapping cis-regulatory chromatin contacts in neural cells links neuropsychiatric disorder risk variants to target genes. Nat Genet. 2019;51:1252–62.
Qin Y, Grimm SA, Roberts JD, Chrysovergis K, Wade PA. Alterations in promoter interaction landscape and transcriptional network underlying metabolic adaptation to diet. Nat Commun. 2020;11:962.
Yue F, Cheng Y, Breschi A, Vierstra J, Wu W, Ryba T, et al. A comparative encyclopedia of DNA elements in the mouse genome. Nature. 2014;515:355–64.
Spicuglia S, Vanhille L. Chromatin signatures of active enhancers. Nucleus. 2012;3:126–31.
Bernstein BE, Mikkelsen TS, Xie X, Kamal M, Huebert DJ, Cuff J, et al. A bivalent chromatin structure marks key developmental genes in embryonic stem cells. Cell. 2006;125:315–26.
Rada-Iglesias A, Bajpai R, Swigut T, Brugmann SA, Flynn RA, Wysocka J. A unique chromatin signature uncovers early developmental enhancers in humans. Nature. 2011;470:279–83.
Cairns J, Orchard WR, Malysheva V, Spivakov M. Chicdiff: a computational pipeline for detecting differential chromosomal interactions in Capture Hi-C data. Bioinformatics. 2019;35:4764–6.
Sarafi MN, Garcia-Zepeda EA, MacLean JA, Charo IF, Luster AD. Murine monocyte chemoattractant protein (MCP)-5: a novel CC chemokine that is a structural and functional homologue of human MCP-1. J Exp Med. 1997;185:99–109.
Xuan W, Qu Q, Zheng B, Xiong S, Fan G-H. The chemotaxis of M1 and M2 macrophages is regulated by different chemokines. J Leukoc Biol. 2015;97:61–9.
Rouault JP, Rimokh R, Tessa C, Paranhos G, Ffrench M, Duret L, et al. BTG1, a member of a new family of antiproliferative genes. EMBO J. 1992;11:1663–70.
Busson M, Carazo A, Seyer P, Grandemange S, Casas F, Pessemesse L, et al. Coactivation of nuclear receptors and myogenic factors induces the major BTG1 influence on muscle differentiation. Oncogene. 2005;24:1698–710.
Stik G, Vidal E, Barrero M, Cuartero S, Vila-Casadesús M, Mendieta-Esteban J, et al. CTCF is dispensable for immune cell transdifferentiation but facilitates an acute inflammatory response. Nat Genet. 2020;52:655–61.
Rodríguez-Ubreva J, Ciudad L, Gómez-Cabrero D, Parra M, Bussmann LH, Di Tullio A, et al. Pre-B cell to macrophage transdifferentiation without significant promoter DNA methylation changes. Nucleic Acids Res. 2012;40:1954–68.
He B, Zhu I, Postnikov Y, Furusawa T, Jenkins L, Nanduri R, et al. Multiple epigenetic factors co-localize with HMGN proteins in A-compartment chromatin. Epigenetics Chromatin. 2022;15:23.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078–9.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923–30.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.
Bunn AG. A dendrochronology program library in R (dplR). Dendrochronologia (Verona). 2008;26:115–24.
Liao Y, Wang J, Jaehnig EJ, Shi Z, Zhang B. WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs. Nucleic Acids Res. 2019;47:W199-205.
Krueger F, Andrews SR. SNPsplit: Allele-specific splitting of alignments between genomes with known SNP genotypes. F1000Res. 2016;5:1479.
Wingett S, Ewels P, Furlan-Magaril M, Nagano T, Schoenfelder S, Fraser P, et al. HiCUP: pipeline for mapping and processing Hi-C data. F1000Res. 2015;4:1310.
Yang T, Zhang F, Yardımcı GG, Song F, Hardison RC, Noble WS, et al. HiCRep: assessing the reproducibility of Hi-C data using a stratum-adjusted correlation coefficient. Genome Res. 2017;27:1939–49.
Durand NC, Robinson JT, Shamim MS, Machol I, Mesirov JP, Lander ES, et al. Juicebox provides a visualization system for Hi-C contact maps with unlimited zoom. Cell Syst. 2016;3:99–101.
Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38:576–89.
Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29:24–6.
Wolff J, Bhardwaj V, Nothjunge S, Richard G, Renschler G, Gilsbach R, et al. Galaxy HiCExplorer: a web server for reproducible Hi-C data analysis, quality control and visualization. Nucleic Acids Res. 2018;46:W11–6.
Cairns J, Freire-Pritchett P, Wingett SW, Várnai C, Dimond A, Plagnol V, et al. CHiCAGO: robust detection of DNA looping interactions in Capture Hi-C data. Genome Biol. 2016;17:127.
Li D, Purushotham D, Harrison JK, Hsu S, Zhuo X, Fan C, et al. WashU epigenome browser update 2022. Nucleic Acids Res. 2022. https://doi.org/10.1093/nar/gkac238.
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.
Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9:R137.
Shao Z, Zhang Y, Yuan G-C, Orkin SH, Waxman DJ. MAnorm: a robust model for quantitative comparison of ChIP-Seq data sets. Genome Biol. 2012;13:R16.
Yates AD, Achuthan P, Akanni W, Allen J, Allen J, Alvarez-Jarreta J, et al. Ensembl 2020. Nucleic Acids Res. 2020;48:D682–8.
Van Oevelen C, Collombet S, Lepoivre C, Bussmann L, Kallin EM, Vicent G, et al. 2015. Genome-wide maps of chromatin states during B-cell to Macrophage lineage switching. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE53173.
Van Oevelen C, Collombet S, Lepoivre C, Bussmann L, Kallin EM, Vicent G, et al. 2015. Genome-wide maps of transcription factor and co-factor binding during B-cell to Macrophage lineage switching. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE53362.
Cellular reprogramming is driven by widespread rewiring of promote-enhancer interactions. NCBI GEO accession GSE221737. 2023. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE221737.
Acknowledgements
We thank the members of the Fraser lab for critical discussions and Thomas Graf for providing the C10 cell line.
Funding
This work was supported by the Frank and Yolande Fowler Foundation.
Author information
Authors and Affiliations
Contributions
BH conducted experiments with the help of YH. MW performed the bioinformatics analysis with help of JS. MW, DS, BH, and PF wrote the manuscript. PF directed the experiments. All the authors read, edited, and approved the manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1:
Fig. S1. Transcriptional profile changes during Pre-B transdifferentiation. a Principal component analysis (PCA) RNAseq of each replicate from different time points. b SNPsplit identifies the percentage of total RNA reads that align to Cebpa gene region that contain Mus musculus specific SNP (mouse specific), Rat specific SNP (Rat specific), regions that do not contain any SNP (Unassignable), conflicting SNP information (Conflicting). c Volcano plots show upregulated and downregulated differentially expressed genes (DEGs). d GO ontology 5 Fuzzy c-means clusters of DEGs e qPCR analysis of Pre-B cell and macrophage marker gene expression during Pre-B cell trans-differentiation. f Gene expression level (TPM) of Pre-B cell-specific genes (Blnk, Cd19, Cd79a, Cd79b, Vpreb1, Vpreb2, Vpreb3), and macrophage-specific genes (Ccl6, Ctsc, Mmp8, Msr1) during Pre-B trans-differentiation. Fig. S2. Hi-C data analysis. aHiCUP pipeline analysis HiC and PCHi-C data, including the following three categories: Cis < 10Kbp, Cis > 10Kbp, and trans-interactions across four-time points. b Heatmaps show Hi-C matrix at 100kb and 25 kb resolution. Fig. S3. Altered TAD and A/B compartment during Pre-B transdifferentiation. a Number of TADs and size distribution at each time point. b GO enrichment analysis of genes changed TAD boundaries from 0h to 48h. c Pearson correlation coefficient of PC1 values (50kb resolution) of the entire genome. d Gene expression level (TPM) in the A and B compartments at four different time points. A two-sided unpaired t-test was performed for the significance test. ePercentage of compartment change for DEGs (differentially expressed genes). f Log2 fold change of gene expression in dynamic compartments. stable (n= 19157), A to B (n= 450), and B to A (n= 474). g Gene ontology analysis of 2105 genes at 48h that switch from A to B and 1552 genes at 48h that switch from B to A. A two-sided unpaired t-test was performed for the significance test. Fig. S4. PCHi-C data analysis. a HiCUP pipeline analysis of PCHi-C data: Cis <10Kbp, Cis >10Kbp, and trans-interactions across four-time points. b Proportion of promoter-promoter (P-P) and promoter-other interactions (P-O) interactions as deduced from PCHi-C. c Expression values (log2(TPM+1)) of genes/promoters classified into 4 groups (at least one PIR overlaps with histone markers). H3K27ac: interaction between gene/promoter and PIRs only with H3K27ac marks. 0h (n=4829); 1h (n=5710);12h (n=5257); 48h (5735). H3K27Ac+H3K27me3: interaction between gene/promoter and PIRs with both H3K27ac and H3K27me3 marks. 0h (n=8075); 1h (n=7044); 12h (n=7926); 48h (n= 7727). H3K27me3: interaction between gene/promoter and PIRs only with H3K27me3. 0h (n=2757);1h (n=2286); 12h (n= 2545); 48h (n=2165). No marks: interaction between gene/promoter and PIRs without any histone modification. 0h (n=3667), 1h (n=4805); 12h (n=4403); 48h (n= 4291). The box plot represents 25 and 75 percentiles with the median. An unpaired two-sided t test is performed for the significance test. d Boxplot shows gene expression values (log2(TPM +1)) of promoters that do not contact any PIRs with H3K27ac categorized into four groups according to the number of PIRs marked by H3K27me3 at four different time points. 0: Promoters without any PIRs marked by H3K27me3. 0h (n=2551), 1h (n=2991); 12h (n=2711); 48h (n=2682). 1: Promoters with only one PIR marked by H3K27me3 0h (n=1247), 1h (n=1100); 12h (n=1188); 48h (n= 955). 1-5: Promoters with one to five PIRs marked by H3K27me3. 0h (n=749), 1h (n=672); 12h (n=697); 48h (n=567). >5: Promoters with more than five PIRs marked by H3K27me3. 0h (n=80), 1h (n=63); 12h (n=103); 48h (n=95). An unpaired two-sided t-test is performed for the significance test. The box plot represents 25 and 75 percentiles with the median. Fig. S5. Integrated analysis of promoter interactions and epigenetic features. a Maf1 and Ssbp1 region that shows promoter interaction changes accompanied by H3K27me3(green), H3K27ac (yellow), and H3K4me1 at interacting regions, along with gene expression TPM of Maf1 and Ssbp1 during Pre-B transdifferentiation. Multiple Maf1 promoters in two HindIII digested fragments b Log2 fold change of gene expression based on the number of H3K4me1 peaks on the increased PIRs (48h vs 0h). The number of histone mark peaks is classified into three groups. <=1: Genes that gain chromatin interactions with less than one histone mark peaks (n=360). 1-5: Genes that gain chromatin interactions with one to five histone mark peaks (n=313). >5: Genes that gain chromatin interactions with greater than five histone mark peaks (n=69). The unpaired two-sided t-test was performed to test significant differences between groups. The box plot represents 25 and 75 percentiles with the median. Fig. S6. Decreased promoter interactions at 48h. a GO term for genes that lost promoter interactions at 48h. b Examples of significantly decreased B cell-specific gene promoter interactions (Pax5,and Clcf1) across time points. Fig. S7. C/EBPα ChIPseq data analysis. a The number of identified C/EBPα ChIP peaks at 1h, 1h, 12 h and 48h. b Comparison of 1h, 12,h or 48h C/EBPα ChIP peaks signal to 0h. M value represents log2 fold change of read density between two samples, and A value represents the average signal intensity. The red dots represent induced C/EBPα binding sites at 1h, 12h, and 48h time points. The blue dots represent C/EBPα binding sites at 0h. The grey dots represent unchanged C/EBPα binding sites at four time points. Fig. S8. Uncropped blots.
Additional file 2:
Table S1.
Additional file 3:
Table S2.
Additional file 4:
Table S3.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Wang, M., He, B., Hao, Y. et al. Cellular reprogramming is driven by widespread rewiring of promoter-enhancer interactions. BMC Biol 21, 264 (2023). https://doi.org/10.1186/s12915-023-01766-0
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12915-023-01766-0