Zixiang Yan, Luzhang Ji, Xiangru Huo, Qianfeng Wang, Yuwen Zhang Bo Wen,2,*

1 MOE Key Laboratory of Metabolism and Molecular Medicine, Department of Biochemistry and Molecular Biology, School of Basic Medical Sciences, and Institutes of Biomedical Sciences, Fudan University, Shanghai 200032, China

2 State Key Laboratory of Genetic Engineering, Collaborative Innovation Center of Genetics and Development, Fudan University,Shanghai 200438, China

KEYWORDS H3K9me2;G9a/GLP;Chromatin organization;3D genome;Genomic compartment

Abstract In the nucleus,chromatin is folded into hierarchical architecture that is tightly linked to various nuclear functions. However, the underlying molecular mechanisms that confer these architectures remain incompletely understood. Here, we investigated the functional roles of H3 lysine 9 dimethylation (H3K9me2), one of the abundant histone modifications, in three-dimensional (3D)genome organization. Unlike in mouse embryonic stem cells, inhibition of methyltransferases G9a and GLP in differentiated cells eliminated H3K9me2 predominantly at A-type(active)genomic compartments, and the level of residual H3K9me2 modifications was strongly associated with B-type (inactive) genomic compartments. Furthermore, chemical inhibition of G9a/GLP in mouse hepatocytes led to decreased chromatin-nuclear lamina interactions mainly at G9a/GLP-sensitive regions,increased degree of genomic compartmentalization,and up-regulation of hundreds of genes that were associated with alterations of the 3D chromatin. Collectively, our data demonstrated essential roles of H3K9me2 in 3D genome organization.

Introduction

The H3 lysine 9 dimethylation (H3K9me2) is one of the histone modifications involved in gene silencing and chromatin repression [1]. Unlike H3K9me3 that associates with constitutively pericentromeric heterochromatin, H3K9me2 marks broad chromatin domains that are named large organized chromatin K9 modifications (LOCKs) across the genome[2,3]. H3K9me2 is mainly deposited by methyltransferases G9a (EHMT2) and GLP (EHMT1), which primarily form a stable heterodimer in vivo. G9a and GLP cannot compensate each other, and depletion of either can significantly reduce the level of H3K9me2; a double knockout, however, does not get further effect[4,5].H3K9 methylation is selectively recognized and bound by heterochromatin protein 1(HP1)family members, HP1α, HP1β, and HP1γ (also known as CBX5,CBX1, and CBX3, respectively). They function to promote heterochromatin formation and maintenance together with other proteins,including G9a[6,7]. Some small-molecule inhibitors, including BIX-01294, UNC0638, and UNC0642, have been developed to specifically inhibit the catalytic activity of G9a/GLP,which provides a powerful tool to study this histone modification [8–10].

It has been well recognized that the heterochromatin is darkly stained and distributed at the nuclear periphery based on the electron microscopy [11]. Through Lamin B1 (LB1)DamID technology, thousands of lamina-associated domains(LADs)have been identified in the mammalian genomes;these LADs dynamically approach nuclear lamina(NL)and contain fewer genes that are primarily repressed and late replicating[12,13]. Interestingly, the genomic regions of large H3K9me2 domains are highly overlapped with those of LADs [2]. Furthermore,via m6A-tracer and fluorescence in situ hybridization(FISH)-based imaging, G9a-mediated H3K9me2 has been proved to be required for chromatin-NL interactions [14,15].

Besides large chromatin domains such as LADs, the genome is folded into multi-scale higher-order architecture in the eukaryotic nuclei. Each chromosome occupies its distinct territory in the nuclear space[16].Based on chromatin conformation capture (3C)-based technologies, the 3D chromatin architecture has been revealed on different levels ranging from 100 kb to 10 Mb, including chromatin loops, topologically associating domains (TADs), and genomic compartments[17–19]. At the megabase scale, genomic compartments are segregated into two types, named A and B, which represent active and inactive chromatin domains,respectively[17].Interestingly, B compartments are correlated with LADs [11,20].

Based on the links among H3K9me2,LADs,and B compartments,H3K9me2 could play roles in 3D genome organization.In hematopoietic progenitors,it has been found that reduction of H3K9me2 by inhibiting G9a/GLP can change chromatin accessibility,as identified with the formaldehyde-assisted isolation of regulatory elements(FAIRE)technique[21].However,details as to whether and how H3K9me2 may alter 3D chromatin architecture remain unknown. Hence, in the present study,we used Hi-C,DamID,ATAC-seq, and ChIP-seq technologies to study the functional roles of H3K9me2 in 3D chromatin organization at the genome-wide scale.

Results

G9a/GLP inhibition removes H3K9me2 mainly at A compartments or inter-LADs in differentiated cells

To investigate the functional roles of H3K9me2 in 3D genome organization,we used mouse AML12 hepatocytes as the model,as our previous studies showed that nuclear architecture was well preserved in AML12 cells [22,23]. To remove H3K9me2,we treated AML12 cells with UNC0638, a selective inhibitor of G9a/GLP[9].We observed a reduction in the H3K9me2 level after UNC0638 treatment by immunofluorescence (IF) and Western blotting(WB)(Figure 1A,Figure S1A and B).It should be noted that the UNC0638 treatment cannot eradicate H3K9me2 completely, and some H3K9me2 signals remained at the nuclear periphery (Figure 1A). We also examined other modifications on the H3K9 residue and observed a slight decrease of H3K9me3 by WB (Figure S1A and B), which was similar to the result previously reported [24]. This observation may be on account of the‘‘binary switch”on the same residue and non-processive H3K9 methylation mechanism by methyltransferases such as SUV39H1 [25,26]. However, the nuclear distribution and genome-wide location of H3K9me3 and H3K9ac (H3 lysine 9 acetylation) were hardly changed, as revealed by IF and ChIP-seq(Figure S1C and D).Additionally,we examined other covalent histone modifications (H3K27ac,H3K27me3, H3K4me1, and H3K4me3) by WB and IF, and found no noticeable effect(Figure S1E and F).

Next, in order to investigate the profile of genome-wide H3K9me2 upon UNC0638 treatment in AML12 cells, we performed quantitative ChIP-seq experiments by spiking in human HeLa cells. The degree of H3K9me2 reduction along the genome was not random as shown by ChIP-seq data,but presented a region-specific pattern (Figure 1B). It should be noted that H3K9me2 ChIP-seq tracks showed more or less equal distribution between LADs and inter-LADs(iLADs) when the data were not normalized. However, after normalizing the ChIP-seq data to input, the average H3K9me2 profile mainly decreased at iLADs or A compartments, but displayed a less degree of reduction at LADs or B compartments (Figure 1C). This result was in agreement with the observation by IF imaging.

To examine H3K9me2 reduction after UNC0638 treatment at higher resolution, we called subcompartments in AML12 cells based on a reported method[19]by integrating Hi-C contact matrices, histone modifications, and lamina-associating maps(Figure S1G and H).The degree of H3K9me2 reduction in A1 and A2 subcompartments were higher than those in B1,B2, and B3 subcompartments (Figure S1I). Moreover,H3K9me2 in B2 subcompartments, which presented a higher H3K9me3 level, were more resistant to UNC0638 treatment than that in B1 and B3 subcompartments (Figure S1I).

Consistently,a reanalysis of published ChIP-seq data of G9a knockout and wild-type mouse embryonic fibroblast (MEF)cells[27]also showed a specific reduction of H3K9me2 at iLADs(Figure S2A).Furthermore,G9a depletion in mouse cardiomyocytes in vivo resulted in similar changes in H3K9me2 patterns[28](Figure S2B).It should be noted that H3K9me2‘‘increased”at LADs in G9a/GLP knockout samples,as these datasets were not generated by using the spike-in method and could not be quantitative. These results indicate that neither chemical nor genetic inhibition of G9a/GLP completely removes H3K9me2 modification at LADs/B compartments.

Considering the different removal patterns of H3K9me2 after G9a/GLP inhibition, we then referred these dramatically H3K9me2-decrease chromatin regions as G9a/GLP-sensitive regions (GSRs) (Figure 1D), which occupied 37.5% of the genomic regions in AML12 cells (Table S1). There are 68.87% and 74.86% of GSRs locating within iLADs and A compartments in AML12 cells,respectively(Figure 1E).These results suggest that the sensitivity of H3K9me2-modified regions upon G9a/GLP inhibition demarcates two different types of genomic compartments, which are analogous to A/B compartments or iLADs/LADs.

G9a/GLP inhibition removes H3K9me2 globally in mESCs

Figure 1 Region-specific removal of H3K9me2 modifications upon G9a/GLP inhibitionA. Representative IF images of H3K9me2 in DMSO and UNC0638 treated AML12 cells. Scale bar, 10 μm. B. Representative ChIP-seq tracks of H3K9me2 upon UNC0638 treatment in AML12 cells as well as PC1 values and smoothed LB1 DamID signals of AML12 cells[22]. Below: zoom-in view of the highlighted region. C. Aligned H3K9me2 profiles of mirrored border regions of LADs (left) and A/B compartments(right)in DMSO and UNC0638 treated AML12 cells.D.Box plots showing relative H3K9me2 levels(DMSO–UNC0638)at GSRs and non-GSRs in AML12 cells. The numbers at the bottom show the genomic bins of 1 kb. E. Percentages of GSRs in LADs/iLADs (upper) and A/B compartments (below) in AML12 cells. IF, immunofluorescence; LB1, Lamin B1; LAD, lamina-associated domain; iLAD, inter-LAD; GSR, G9a/GLP-sensitive region.

As H3K9me2 modification is highly dynamic during stem cell differentiation [2,10], we wondered whether the sensitivity of G9a/GLP inhibition on H3K9me2 was cell-type specific.Treating mouse embryonic stem cells (mESCs) with UNC0638 induced decrease of H3K9me2 (Figure S3A).Unlike the results observed in AML12 cells, H3K9me2 modifications at both the nuclear interior and periphery were almost removed completely, and there was a foci-like staining of H3K9me2 in mESCs treated with UNC0638 as shown by IF images,which was similar to H3K9me3 staining(Figure S3B).A similar result was also found in G9a-/-mESCs[29].Consistently, ChIP-seq data also showed that most of H3K9me2 modifications,including those in LADs/B compartments,were removed entirely, but some ‘‘sharp peaks” remained (Figure S3C and D). These retained H3K9me2 peaks were mainly(70.3%) overlapped with H3K9me3 peaks (Figure S3E). We also analyzed the features of the remaining non-overlapping peaks, and found that they were mostly intergenic and enriched with repeat sequences (Figure S3F and G). In addition, the ChIP-seq data of GLP knockout mESCs [27] yielded similar results(Figure S3C).Therefore,unlike in differentiated cells such as AML12 and MEFs, H3K9me2 modification in mESCs was mainly G9a/GLP-sensitive. Although H3K9me2 was globally removed after G9a/GLP inhibition in mESCs,the A/B compartments were largely maintained (Figure S3H).We thus speculated that those remaining peaks represented the more condensed heterochromatin, which is regulated by other methyltransferases independent of G9a/GLP.

H3K9me2 insensitive to G9a/GLP inhibition is associated with inactive genomic compartments

As shown by the representative region in Figure 2A,the tracks of LB1 DamID and PC1 displayed similar wavy patterns in AML12 cells; after UNC0638 treatment, H3K9me2 was more reduced at regions with weaker LB1 DamID signals and stronger A-tendency PC1 values, in both iLADs/A and LADs/B;the wavy patterns of remaining H3K9me2 were more similar to those of LB1 DamID and PC1 tracks.Quantitative analysis showed that, the H3K9me2 levels after G9a/GLP inhibition exhibited much stronger correlation with PC1 values(r = -0.661) or LB1 DamID signals (r = 0.701), comparing with the DMSO-treated samples (r = -0.198 for PC1,r=0.205 for LB1 DamID;Figure 2B and C).However,other inactive histone modifications, such as H3K9me3 and H3K27me3, displayed little correlation with PC1 or LB1 DamID (Figure 2D and E).

Figure 2 G9a/GLP-insensitive H3K9me2 is associated with compartmental statusA. Representative tracks of H3K9me2 ChIP-seq and difference (DMSO – UNC0638) in DMSO and UNC0638 treated AML12 cells as well as PC1 values and smoothed LB1 DamID signals of AML12 cells[22].Below:zoom-in view of the highlighted region.B.Correlation between H3K9me2 levels and PC1 values in DMSO (left)and UNC0638 (right)treated AML12 cells. C.Correlation between H3K9me2 levels and LB1 DamID signals in DMSO (left) and UNC0638 (right) treated AML12 cells. D. Correlation between PC1 values and H3K9me3 (left)/H3K27me3 (right) levels in DMSO treated AML12 cells. E. Correlation between LB1 DamID signals and H3K9me3(left)/H3K27me3(right)levels in DMSO treated AML12 cells.The 40 kb bins were used for the correlation analyses shown from B to E.

Histone modifications can reflect different chromatin state segments, which are associated with fine-scale genomic compartmentalization [30]. Our data further suggest that the intrinsic H3K9me2 modifications at non-GSRs can better reflect the inactive environment of chromatin and are highly correlated with the genomic compartmentalization.

G9a/GLP inhibition decreases chromatin-NL interactions primarily at GSRs

Figure 3 Inhibition of G9a/GLP by UNC0638 treatment decreases chromatin-NL interactions mainly at GSRsA. Representative electron micrographs of DMSO and UNC0638 treated AML12 cells. Black arrowheads indicate the nuclear periphery chromatin.Scale bar,2 μm.B.Representative tracks of smoothed LB1 DamID signals and H3K9me2 ChIP-seq data of DMSO and UNC0638 treated AML12 cells.Black bars represent locations of LADs,and colored shadows indicate switched regions.C.Proportions of switched and stable LADs inAML12 cellsupon UNC0638 treatment.D.Box plot showingrelative expressionchangesof genesassociated with loss orgain of LADs in AML12 cells upon UNC0638 treatment.P values,Wilcoxon rank-sum test.E.Proportions of the four sections(GSRs/LADs,GSRs/iLADs,non-GSRs/LADs,and non-GSRs/iLADs)in the genome of AML12 cells.F.Box pots of LB1 DamID signals of the four sections in DMSO and UNC0638 treated AML12 cells.The numbers above or below the box plots show the median of plots.The numbers at the bottom show the probe counts of LB1 DamID.P values,Student’s t-test.NL,nuclear lamina;CP,cytoplasm;NP,nucleoplasm.

Our electron microscopy data showed that less nuclear periphery heterochromatin was seen in UNC0638 treated AML12 cells(Figure 3A).To investigate the effects of G9a/GLP inhibition on chromatin-NL interactions at a genome-wide scale,we conducted LB1 DamID assay in AML12 cells. By examining the LB1 DamID signals, the loss or gain of LADs occurred in some genomic regions (Figure 3B), and the coverage of LADs decreased from 54.75%to 47.43%(Figure 3C;Table S2).These switched LADs or iLADs were associated with the significant changes of gene expression levels(Figure 3D).

Figure 4 Removal of H3K9me2 at GSRs increases genomic compartmentalizationA. Representative Hi-C contact matrices (Obs/Exp) in DMSO (left) and UNC0638 (right) treated AML12 cells. B. Average contact enrichment between pairs of 250 kb bins ranked by PC1 values in DMSO (left) and UNC0638 (middle) treated AML12 cells, and the difference between them (right). C. Box plot showing compartmentalization strength across chromosomes in DMSO and UNC0638 treated AML12 cells. The numbers below the boxes show the median. P value, two-factor ANOVA. D. Violin plots showing TAD boundary strength changes within A(left)and B(middle)compartments and between A and B(right)compartments in AML12 cells upon UNC0638 treatment. E. Box plots showing interaction scores of the four fractions (GSRs/LADs, GSRs/iLADs, non-GSRs/LADs, and non-GSRs/iLADs) in DMSO and UNC0638 treated AML12 cells. The numbers above the boxes show the median. P values, two-factor ANOVA. TAD, topologically associating domain.

In order to further study the connection between H3K9me2 and chromatin-NL interactions, we classified the genome into four sections: GSRs/LADs (11.66%), GSRs/iLADs (25.80%),non-GSRs/LADs (41.33%), and non-GSRs/iLADs (21.81%)(Figure 3E). Because of the large sample size in the statistical tests, chromatin-NL interactions were significantly different in all these sections, but the degree and direction of changes were variable.As shown in Figure 3F,the chromatin-NL interactions in GSRs/iLADs and GSRs/LADs decreased after UNC0638 treatment, consistent with the obvious decrease of H3K9me2 in these regions. However, in non-GSRs, where H3K9me2 were less reduced, the chromatin-NL interactions increased (non-GSRs/LADs) or were hardly changed(non-GSRs/iLADs). Of note, the detected increase of LB1 DamID signals in non-GSRs/LADs could be not biologically relevant, but resulted from the normalization of LB1 DamID data of UNC0638 treated samples. Therefore, G9a/GLP inhibition decreases chromatin-NL interactions mainly at GSRs.

G9a/GLP inhibition increases genomic compartmentalization

Figure 5 Gene expression changes upon UNC0638 treatment in AML12 cellsA.Heatmap of differentially expressed genes between DMSO and UNC0638 treated AML12 cells.Two biological replicates(independent DMSO/UNC0638 treatment,RNA library preparation,and sequencing)were performed for each sample.B.Top ten Gene Ontology terms of up-regulated genes upon UNC0638 treatment in AML12 cells.C.Quantification of expression changed/unchanged genes in four sections(GSRs/LADs, non-GSRs/LADs, GSRs/iLADs, and non-GSRs/iLADs) in AML12 cells. The percentages of changed/unchanged genes in each fraction were shown. D.Clusters of gene-centric signal changes of H3K9me2 ChIP-seq signals, PC1 values,LB1 DamID signals and ATAC-seq signals upon UNC0638 treatment in AML12 cells.Relative signal changes were calculated at promoters with bin size of 5 kb.E.Box plot showing gene expression changes in four clusters upon UNC0638 treatment in AML12 cells.P values,Wilcoxon rank-sum test.F.Proportions of the cluster 4 genes in four sections(GSRs/LADs,GSRs/iLADs,non-GSRs/LADs,and non-GSRs/iLADs)in AML12 cells.

To further investigate the impacts of H3K9me2 on 3D genome organization, we generated high-resolution chromatin interaction maps by in situ Hi-C assay in DMSO and UNC0638 treated AML12 cells. We performed two biological replicates for the Hi-C experiments and obtained highly consistent results(Figure S4A–C). Then, by overviewing the long-range Hi-C contact maps, we could see interaction changes among compartments (Figure 4A). Further analysis of genomic compartmentalization by pairwise comparison of chromatin interactions between ranked PC1 intervals revealed that the interactions between A and B compartments decreased, while the interactions between the same type of compartments (A vs. A; B vs. B) increased (Figure 4B, Figure S4D). Through quantification with the duplicates of Hi-C data, we showed that compartmentalization strength [31], increased by 23%,from 3.495 to 4.31, upon G9a/GLP inhibition (Figure 4C).Similarly, the strength of TAD boundaries between A and B compartments significantly increased, but those within A or B compartments were not obviously changed, indicating the increased insulation between A and B compartments (Figure 4D; Table S3). These results suggest the increased degree of genomic compartmentalization upon G9a/GLP inhibition,which is coincident with the enlarged difference of H3K9me2 levels between A and B compartments.

Furthermore,via the stratification of the genome by considering GSRs and LADs, we found that interaction scores,which measure the chromatin interactions at given regions[23], increased significantly in all sections except for non-GSRs/iLADs after UNC0638 treatment (Figure 4E). These observations were consistent with the results of chromatin-NL interactions and in agreement with the previous report that chromatin-NL was coupled with chromatin-chromatin interactions in the genome organization [32]. In non-GSRs/LADs where H3K9me2 was less removed, the increased chromatinchromatin interactions could be the secondary effect of chromatin changes at GSRs after H3K9me2 reduction. Similarly,it has been reported that active chromatin marks on euchromatin can drive spatial sequestration of heterochromatin indirectly in differentiated cells [33].

H3K9me2 reduction by G9a/GLP inhibition alters the expression of hundreds of genes

To further investigate the biological functions of G9a/GLPsensitive H3K9me2 in AML12 cells, we treated the cells with UNC0638 and examined gene expression by RNA-seq. There were 710 genes with significant expression changes upon UNC0638 treatment (Table S4). Among them, 484 genes(68%) were up-regulated, and 226 genes (32%) were downregulated (Figure 5A). Gene Ontology (GO) analysis showed that up-regulated genes were mainly enriched in liver biological processes, including steroid biosynthesis and lipid metabolism (Figure 5B), yet down-regulated genes had no significant enrichment.

To explore the relationship between gene expression and chromatin organization, we examined genes in the context of GSRs and LADs. In all the four sections (GSRs/LADs,GSRs/iLADs, non-GSRs/LADs, and non-GSRs/iLADs), the numbers of up-regulated genes were more than those of down-regulated genes (Figure 5C). Usually, gene upregulation is associated with the increased chromatin accessibility. Thus, we performed ATAC-seq to examine the changes of chromatin accessibility upon UNC0638 treatment, and combined its results with the high-throughput data, including H3K9me2 ChIP-seq, RNA-seq, LB1 DamID, and Hi-C, for an integrated analysis. After excluding the genes that were not expressed in both DMSO and UNC0638 treated samples,we used a gene-centric approach to aggregate relative signal changes of these chromatin features and identified four clusters showing different chromatin alterations(Figure 5D;Table S5).In all the four clusters, H3K9me2 more or less decreased, and the ATAC-seq signals increased averagely, indicating the increased chromatin accessibility at promoters of these genes.Among them, cluster 4 showed decreased H3K9me2 level,increased Hi-C PC1(compartment A tendency),and decreased chromatin-NL interactions (LB1 DamID) (Figure 5D). Genes in cluster 4 tended to be up-regulated after UNC0638 treatment; however, the gene expression levels of other three clusters were not obviously affected (Figure 5E). Moreover,genes in cluster 4 were mainly enriched in GSRs (64.89% in GSRs/iLADs and 9.48% in GSRs/LADs; Figure 5F). These results indicate that the removal of G9a/GLP-sensitive H3K9me2 up-regulates the expression of hundreds of genes,which are associated with the more open chromatin environment.

Discussion

In this study, we showed that G9a/GLP inhibition mainly removed H3K9me2 at A compartments. G9a/GLPinsensitive H3K9me2 was highly correlative with inactive genomic compartments globally, suggesting that the level of G9a/GLP-sensitivity of H3K9me2 could reflect the intrinsic tendency of chromosomes to segregate into diverse compartments. Chromatin compartmentalization can be regulated by histone modifications through phase separation [34]. In this study, we showed that the removal of G9a/GLP-sensitive H3K9me2 in AML12 cells could increase genomic compartmentalization and TAD boundary strength between A and B compartments, and decrease chromatin-NL interactions heterogeneously that are associated with H3K9me2 reduction. Additionally, H3K9me2 reduction by G9a/GLP inhibition up-regulated hundreds of genes associated with alterations of chromatin organization. Our results thus revealed the functional roles of H3K9me2 in 3D chromatin organization.

It has been shown that, based on immunostaining data,G9a is mainly responsible for H3K9 methylation at euchromatic loci (i.e., nuclear areas other than the DAPI-dense pericentromeric heterochromatin)in mESCs[29].In this study,we further demonstrated that regardless of chemical inhibition or genetic depletion of G9a/GLP, H3K9me2 was not removed completely in B compartments or LADs in differentiated cells.However, H3K9me2 can be removed more effectively in mESCs by G9a/GLP inhibition. Moreover, the induction of mESCs to epiblast-like cells and primordial germ cell-like cells is accompanied by a large-scale reorganization of chromatin signatures, including H3K9me2 patterns [35]. During hematopoietic stem cell lineage commitment, G9a/GLPdependent H3K9me2 marks are established gradually [10].These data imply a complexity of mechanisms underlying H3K9me2 modification in differentiated cells. The incomplete removal of H3K9me2 by G9a/GLP inhibition in differentiated cells could imply an inefficient UNC0638 treatment at B compartments or a G9a/GLP-independent H3K9me2 modifying system existing in differentiated cells. A recent study has reported that histone turnover of heterochromatin at nuclear periphery is repressed,which may also explain the insensitivity of H3K9me2 to G9a/GLP inhibition in B compartments [36].

However, it remains unclear what chromatin modifiers could be responsible for the remaining H3K9me2 modification at inactive compartments of differentiated cells after G9a/GLP inhibition. The Suv39 family members, including SUV39H1/2 and SETDB1/2, are also known to be involved in H3K9me2 modification [1]. It has been reported that during somatic cell reprogramming into induced pluripotent stem cells (iPSCs),H3K9 methylation is a crucial barrier and knockdown of Suv39h1/2 and Setdb1 can promote the reprogramming, but the H3K9me2 level is not significantly decreased at some gene loci by ChIP-qPCR examination[37].A recent report has also shown that even triple-knockout of Suv39h1/2 and Setdb1 cannot cause a marked decrease of H3K9me2 in mouse liver [38].It would be fascinating to further investigate the modifiers of these G9a/GLP-insensitive H3K9me2 marks and their functions in chromatin organization and other biological processes in the future.

The eukaryotic NL is mainly composed of Lamin A/C and Lamin B, as well as other interacting components inside or nearby the networks,which provides a depressing environment for LADs. G9a is found to interact with the NL-interacting protein BAF [39].Also, G9a and GLP play novel roles in regulating heterochromatin anchorage to the nuclear periphery via the methylation and stabilization of Lamin B1,which associates with H3K9me2-marked peripheral heterochromatin[40].H3K9me3 in LADs has been reported to be combined with Lamin B receptor mediated by CBX5[41].In mammalian cells,the histone deacetylase,HDAC3,has been reported to interact with the Lamin-associated protein LAP2β to maintain the repressed state of peripheral chromatin [42]. In the current study,we provided genomic data showing that G9a/GLP inhibition mainly decreased chromatin-NL interactions of LADs at GSRs,whereas LADs at non-GSRs remained intact.Therefore, for a better understanding of the organization and functions of LADs not affected by G9a/GLP inhibition, further efforts in the field would be needed to decipher factors associated with G9a/GLP-insensitive H3K9me2 at the nuclear periphery.

Materials and methods

Cell culture and UNC0638 treatment

The mouse hepatocyte cell line alpha mouse liver 12(AML12;catalog No. CRL-2254, ATCC, Manassas, VA) was cultured in 37°C and 5%CO2incubator in medium DMEM/F12(Catalog No. 11320033, ThermoFisher Scientific, Waltham, MA)supplemented with 10% fetal bovine serum (FBS; Catalog No. 16140071, Gibco, Grand Island, NY), 1 × ITS Liquid Media Supplement (100×; catalog No. 41400045, Gibco) and 40 ng/ml dexamethasone (Catalog No. D4902, Sigma, Darmstadt, Germany). AML12 cells were treated with 8 μM of UNC0638 (Catalog No. S8071, Selleck, Houston, TX) for 5 days, and the same amount of DMSO (Catalog No. D2650,Sigma) was used as a control.

Mouse embryonic stem cell line E14 (E14TG2a) was cultured in 2i/LIF conditions as described [12]. The G9a/GLP inhibitor UNC0638 was added into the medium to a final concentration of 0.5 μM for 5 days of treatment, and the same amount of DMSO was used as a control.

Western blotting

Cell lysis was incubated at 95 °C for 15 min with 1× SDSPAGE Loading Buffer.The primary antibodies are as follows:anti-H3 (Catalog No.ab1791,Abcam, Cambridge,UK), anti-H3K9me2 (Catalog No. A2359, Abclonal, Woburn, MA),anti-H3K9me3 (Catalog No. A2360, Abclonal), anti-H3K9ac(Catalog No. 39137, Active Motif, Carlsbad, CA), anti-H3K27ac (Catalog No. 39133, Active Motif), anti-H3K27me3 (Catalog No. 9733, Cell Signaling Technology,Boston, MA), anti-H3K4me1 (Catalog No. 39297, Active Motif), and anti-H3K4me3 (Catalog No. A2357, Abclonal).HRP-conjunction secondary antibodies (Catalog No.711-035-152, Jackson ImmunoResearch Laboratories, West Grove, PA) and Chemiluminescent HRP Substrate (Catalog No. WBKLS0500, Millipore, Billerica, MA) were used for the detection. Software AlphaView was used for the relative grayscale statistics.

Immunofluorescence assay

Cells were fixed in 4%formaldehyde followed by the treatment with 0.5%Triton X-100 for 10 min at room temperature.The cells were blocked with phosphate buffered saline (PBS) containing 4% bovine serum albumin (BSA) for 30 min at room temperature and then processed for immunostaining. The primary antibodies are the same as that we used for WB and the secondary antibodies are as follows: Alexa Fluor 488-AffiniPure Donkey Anti-Rabbit IgG (Catalog No.711-545-152, Jackson ImmunoResearch Laboratories) and Alexa Fluor 594-AffiniPure Donkey Anti-Rabbit IgG(Catalog No. 711-585-152, Jackson ImmunoResearch Laboratories).The slides were counterstained with DAPI (Catalog No.C1002,Beyotime,Shanghai,China).Fluorescence images were taken with Leica confocal microscope (Catalog No. TCS SP5,Leica, Wetzlar, Germany), and we controlled the confocal parameters unchanged in each set of experiments. The images were analyzed by the LAS AF Lite software.

Transmission electron microscopy

Cells were scratched and collected followed by fixing in 2.5%glutaraldehyde diluted in 0.1 M phosphate buffer (pH 7.4).The samples were further processed with 1% osmium tetroxide, dehydrated with acetone, and embedded with resin, at the EM facility in the School of Basic Medical Sciences,Fudan University, China. Target cells were randomly selected and captured with a transmission electron microscope (Catalog No. JEM-1230, JEOL, Tokyo, Japan).

ChIP-seq

The ChIP experiments were performed as described [12] with the antibodies anti-H3K9me2 (Catalog No. ab1220, Abcam),anti-H3K9me3(Catalog No.ab8898,Abcam),H3K9ac(Catalog No. 39137, Active motif), and anti-CTCF (Catalog No.3418, Cell Signaling Technology). To calibrate H3K9me2 ChIP-seq, HeLa cells were mixed with DMSO or UNC0638 treated AML12 cells in 1:4 proportion initially. The libraries were prepared using the NEBNext Ultra II DNA Library Prep Kit (Catalog No. E7645, NEB, Ipswich, MA) followed by next-generation sequencing (NGS) using the Illumina HiSeq X Ten system. Two biological replicates (independent DMSO/UNC0638 treatment, ChIP assay, and sequencing)were performed. Raw sequencing reads were mapped to the mouse mm9 genome reference using Bowtie2 [43]. Duplicated read pairs were discarded using Samtools [44]. Log2ratio of IP and input signals was calculated using deepTools [45].For the normalization of spike-in H3K9me2 ChIP-seq data,we referred to the reported method [46]. The NGS data were mapped to refseq genome of mm9(AML12)and hg19(HeLa),respectively. Then, the Input, DMSO treated and UNC0638 treated AML12 H3K9me2 ChIP-seq signals were divided by the ratio of unique mapping reads[human/(human+mouse)].

To define GSRs, we firstly subtracted DMSO and UNC0638 treated ChIP signals with bins of 1 kb. Then, the data were normalized and smoothed using a moving average approach with 40 kb window size. At each bin, we converted the smoothed signals into the t-statistics and then determined a threshold -8 for calling domains based on probe level FDR < 0.01 using the left part of the distribution as null.GSRs were defined as regions with consecutive bins with tstatistics greater than the threshold,and the length of domains was limited to no less than 50 kb. Meanwhile, domains with distances less than 10 kb were merged.

RNA-seq

RNA was extracted with TRIzol Reagent (Catalog No.15596018, Life Technologies, Grand Island, NY). The RNA libraries were prepared using the Ribo-Zero Gold rRNA Removal Reagent (Catalog No. MRZG12324, Illumina, Hayward, CA) and NEBNext Ultra II RNA Directional Library Prep Kit (Catalog No. E7760, NEB) followed by NGS using the Illumina HiSeq X Ten system. Two biological replicates(independent DMSO/UNC0638 treatment and sequencing)were performed. Paired-end sequencing reads were mapped to the mouse mm10 genome reference by HISAT2 [47]. To identify differentially expressed genes, reads mapped to annotated genes were counted with the HTSeq package [48]. Fold change and P values were calculated by DESeq2 [49].FDR<0.05 and|Log2fold change|>1 were used as thresholds. GO term enrichment analysis was performed using the DAVID tool [50].

To define gene-centric changes of chromatin features, we removed genes that were not expressed in both DMSO and UNC0638 treated samples (gene baseMean less than 10.30 in DESeq2) and calculated relative signal changes at promoters with bin size of 5 kb.We used difference(UNC0638-DMSO)for PC1 values and LB1 DamID signals, and ratio (Log2UNC0638/DMSO) for signal changes of ATAC-seq and ChIP-seq. Then we used K-means clustering to divide genes into four categories with different chromatin feature changes.

LB1 DamID

LB1 DamID experiments were performed as described [12].Microarray hybridization assay and data analysis were conducted as described [51]. Two biological replicates (independent DMSO/UNC0638 treatment, LB1 DamID, and microarray hybridization assay) were performed.

ATAC-seq

ATAC experiments and data analysis were conducted as described [52] with minor modifications: harvesting 30,000 AML12 cells followed by transposition reaction at 37 °C for 45 min. Libraries were sequenced via the Illumina HiSeq X Ten system. Two biological replicates (independent DMSO/UNC0638 treatment, ATAC assay, and sequencing) were performed. Paired-end reads were mapped to the mouse genome(mm9) using Bowtie2 with parameters ‘‘-X 2000”. Mitochondrial and PCR duplicate reads were discarded after alignment.MACS2 was used for peak calling[53]with a q value threshold of 0.01.

In situ Hi-C

The in situ Hi-C libraries were prepared as previously described [19]. The libraries were then sequenced via the Illumina HiSeq X Ten system.Two biological replicates were performed (independent DMSO/UNC0638 treatment, Hi-C assay, and sequencing). Paired-end reads were aligned to mm9 reference genome using HiC-Pro [54].

To call compartments,eigenvector function from juicer[55]was used to calculate PC1 values at 40 kb resolution.The sign of PC1 values was adjusted based on gene density to assign A and B compartments.

To call TADs, HiCtool [56] based on the directionality index and hidden Markov model was used at 40 kb resolution.To compare the structure of TADs between two conditions,‘‘the distance between the centers of two boundaries less than or equal to 40 kb” was used as a criterion. To compare TAD boundary strength, the script matrix2insulation.pl [57] was used from cworld::dekker.

To compare interactions between two conditions, the ratio of observed and expected normalized counts was calculated by dump function from juicer.

To calculate compartmentalization strength, we used median (AA, BB)/median (AB) as a measure, where AA is Obs/Exp between pairs of loci with a strong A compartment signal(top 20%based on PC1 values),BB is Obs/Exp between pairs of loci with a strong B compartment signal(bottom 20%based on PC1 values),and AB is Obs/Exp between pairs of loci with the strong A and B compartment signals.

To calculate interaction scores,we analyzed them at the bin level of 40 kb resolution according to their PC1 values. For each bin, we calculated its average interaction frequency with bins of the same section (Cx, the same section means GSRs/LADs, GSRs/iLADs, non-GSRs/LADs, or non-GSRs/iLADs) and its average interaction with any other bins(Ctotal) using Obs/Exp. We calculated the interaction score as CSx= Cx/Ctotal.

Data availability

All the NGS data of LB1 DamID, ATAC-seq, RNA-seq,ChIP-seq, and Hi-C from this study have been deposited in the Genome Sequence Archive [58] at the National Genomics Data Center,Beijing Institute of Genomics,Chinese Academy of Sciences/China National Center for Bioinformation(GSA:CRA002762),and are publicly accessible at https://bigd.big.ac.cn/gsa.

CRediT author statement

Zixiang Yan: Conceptualization, Methodology, Validation,Formal analysis, Investigation, Data curation, Visualization,Writing - original draft. Luzhang Ji: Methodology, Software,Validation, Formal analysis, Data curation, Visualization,Writing - original draft. Xiangru Huo: Validation, Investigation, Writing - original draft. Qianfeng Wang: Investigation.Yuwen Zhang: Investigation. Bo Wen: Conceptualization,Resources, Writing - review & editing, Supervision, Project administration, Funding acquisition. All authors read and approved the final manuscript.

Competing interests

The authors declared that they have no conflict of interest.

Acknowledgments

We thank members of the Wen lab for their suggestions and supports,Jin Li and Yalin Huang for assistance with confocal microscope operation, and Shuangshuang Dong from Feng Zhang Lab for assistance with microarray assay of DamID.This work was supported by the National Key R&D Program of China (Grant No. 2018YFC1003500 awarded to BW) and the National Natural Science Foundation of China (Grant No. 31771435 awarded to BW).

Supplementary material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.gpb.2020.08.001.

ORCID

0000-0002-3278-2479 (Zixiang Yan)

0000-0002-3345-5533 (Luzhang Ji)

0000-0002-7466-4665 (Xiangru Huo)

0000-0001-5107-1595 (Qianfeng Wang)

0000-0002-1227-8427 (Yuwen Zhang)

0000-0003-0394-7366 (Bo Wen)