Asymmetric distribution of parental H3K9me3 in S phase silences L1 elements

Abstract

In eukaryotes, repetitive DNA sequences are transcriptionally silenced through histone H3 lysine 9 trimethylation (H3K9me3). Loss of silencing of the repeat elements leads to genome instability and human diseases, including cancer and ageing^1-3. Although the role of H3K9me3 in the establishment and maintenance of heterochromatin silencing has been extensively studied^4-6, the pattern and mechanism that underlie the partitioning of parental H3K9me3 at replicating DNA strands are unknown. Here we report that H3K9me3 is preferentially transferred onto the leading strands of replication forks, which occurs predominantly at long interspersed nuclear element (LINE) retrotransposons (also known as LINE-1s or L1s) that are theoretically transcribed in the head-on direction with replication fork movement. Mechanistically, the human silencing hub (HUSH) complex interacts with the leading-strand DNA polymerase Pol ε and contributes to the asymmetric segregation of H3K9me3. Cells deficient in Pol ε subunits (POLE3 and POLE4) or the HUSH complex (MPP8 and TASOR) show compromised H3K9me3 asymmetry and increased LINE expression. Similar results were obtained in cells expressing a MPP8 mutant defective in H3K9me3 binding and in TASOR mutants with reduced interactions with Pol ε. These results reveal an unexpected mechanism whereby the HUSH complex functions with Pol ε to promote asymmetric H3K9me3 distribution at head-on LINEs to suppress their expression in S phase.

Extended Data Figure 1.. Leading strand bias of H3K9me3 eSPAN is detected by H3K9me3 antibodies from three different sources.

a. The TA skew and average BrdU-IP-ssSeq bias (BrdU bias) around replication origins in mES cells. b. The TA skew was correlated with BrdU bias. The BrdU bias, reflecting the relative amount of DNA synthesis at leading and lagging strands, was calculated using formula $(\mathrm{W}-\mathrm{C})∕(\mathrm{W}+\mathrm{C})$. $\mathrm{W}$ and $\mathrm{C}$ represent sequencing reads of Watson and Crick strands, respectively. Spearman’s rank correlation coefficient was shown. Each dot represents a 1 kb bin within the 1,928 initiation zones in mES cells. p < 2.2e-16. c. Average eSPAN bias of MCM2, a subunit of the CMG replicative helicase, and 7 histone modifications (H4K20me2, H3K36me3, H3K36me2, H3K4me3, H4K12ac, H4K5ac and H4K5/K12ac) in mES cells. Two independent repeats, indicating by blue and red lines, for each eSPAN experiment shown. d. Raw average H3K9me3 eSPAN bias (ctr) and after normalizing against BrdU bias (no BrdU bias) or TA skew (no TA skew) in mES cells. e. Average bias of H3K9me3 eSPAN bias generated using three H3K9me3 antibodies from different sources. f. Immunoblots of H3K9me3 in different amounts of mES cell lysates by three different H3K9me3 antibodies used in e. Recombinant histone H3/H4 (re.) were used as negative controls. * indicates non-specific signals detected by antibody 2 (Ab2, Diagenode) and antibody 3 (Ab3, Active motif) after heavy exposure. Ab1: self-made in the laboratory and used in this study. n = 3. g. Genome-wide correlations of ENCODE H3K9me3 ChIP-seq with H3K9me3 CUT&Tag signals generated with three different antibodies Ab1, Ab2 and Ab3, with a window size of 20 kb. Note that Ab1 showed the strongest correlation with published ChIP-seq datasets, consistent with the immunoblotting results. These differences in performances of the three H3K9me3 antibodies likely contribute to the different H3K9me3 eSPAN biases observed in e. h. A snapshot of ENCODE H3K9me3 ChIP-seq signals and two repeats of H3K9me3 CUT&Tag signals generated by Ab1 H3K9me3 antibodies at the indicated mouse Chr14 region. For gel source data, see Supplementary Figure 1.

Extended Data Figure 2.. The enrichment of H3K9me3 at the leading strands is also detected in HeLa and primary mouse B cells.

a. Normalized density of ATAC-seq signals around replication origins in mES cells. Black and red lines indicate two independent datasets. b. Raw average H3K9me3 eSPAN bias (ctr) and after removing eSPAN sequencing reads at regions that also contain ATAC-seq peaks (no ATAC) from analysis. c. Correlations between H3K9me3 eSPAN bias and normalized published ATAC-seq signals. Each dot represents a 1 kb bin within the 1,928 initiation zones in mES cells. Spearman’s rank correlation coefficient and p value were shown. d. OK-seq biases at origins in mES (n = 1,928), HeLa (n = 2,809) and primary mouse B cells (n = 1,073) used in this study. e, f. Heatmaps of eSPAN biases of H3K9me3, H3K27me3 and H4K20me2 and OK-seq bias in HeLa (e) and activated mouse B cells (f) at each individual replication origin, with the number of origins used for analysis shown. The heatmap was sorted based on replication efficiency defined by OK-seq. g, h. A snapshot of ChIP-seq and CUT&Tag signals and calculated eSPAN bias for H3K9me3 and H3K27me3 in HeLa (g) and activated mouse B cells (h). OK-seq bias indicates origin location and DNA replication direction (shown by arrow), and L1 elements (≥ 1 kb) with their transcription direction at each locus were shown.

Extended Data Figure 3.. LINE retrotransposons contribute to the asymmetric H3K9me3 distribution.

a. Correlations between OK-seq bias and eSPAN bias of H3K9me3, H4K20me2, H3K36me3, or H3K27me3. Spearman’s rank correlation coefficient and the density distribution were shown. b. Experimental schemes for eSPAN analysis in synchronized mES cells shown in Figure 2a and Extended Data Figure 3c. After pulsing cells with BrdU for 30 min, cells were either sorted by flow cytometry based on the FUCCI reporters (Figure 2a) or treated with nocodazole (Extended Data Figure 3c). See Materials and Methods section for more details. c. Average H3K9me3 eSPAN in asynchronized mES cells or cells synchronized at G2/M phase. Bottom: flow cytometry analysis of cell cycle of asynchronized mES cells and cells arrested at G2/M by nocodazole. d. The relative enrichment of different repetitive elements at 1kb bins with high or low H3K9me3 eSPAN bias surrounding DNA replication origins in HeLa cells. DNA sequences around replication origins were fragmented into 1 kb bin, and ranked based on H3K9me3 eSPAN bias. The top 25% of regions with the highest H3K9me3 eSPAN bias and the bottom 25% with the lowest H3K9me3 eSPAN bias were then used for calculating the enrichment of each indicated DNA element. Fold enrichment is defined as the ratio between the calculated and expected enrichment. com., complexity; rep., repeats. e. Percentage of the accumulative H3K9me3 ChIP-seq signals at different TEs around replication origins with highest (top quartile) and lowest (bottom quartile) H3K9me3 eSPAN bias defined in Figure 2b. f. A schematic representation showing L1 elements whose transcription direction is head-on (HO) and co-direction (CD) with the direction of replication fork movement (left). The numbers and average H3K9me3 eSPAN bias of different TEs within the (−100 kb to 100 kb) regions of 1,928 origins in mES cells and 2,809 origins in Hela cells were counted and shown. Others: all other TEs excluding LINEs. g. Box plots of H3K9me3 eSPAN bias at HO L1s separated by their locations in the early, mid or late replicating origins, which are defined based on the replication timing data in mES cells. h. Box plots of H3K9me3 eSPAN bias at HO L1s separated by their locations in genome compartment A or B based on Hi-C datasets in mES cells. Box plots (g, h) show the median, 25% and 75% quartiles and minimal and maximal values with p values by two-sided Mann–Whitney–Wilcoxon tests. Each panel is a representative of at least two independent experiments. See Materials and Methods for more details.

Extended Data Figure 4.. Strong asymmetric H3K9me3 distribution is detected at “young” and long L1s.

a. All HO L1 families with more than 195 copies in mES (left) and HeLa cells (right) were ranked from left to right based on average H3K9me3 eSPAN bias. See Source Data for numbers of L1s in each family. Data were plotted as mean ± SD. Note that L1Md_T and L1Md_A in mES and L1PA in HeLa cells, which have been reported as young and full-length L1s in mouse and human, respectively, show bigger bias. b. Pearson correlations between H3K9me3 eSPAN bias at HO L1s surrounding origins and their corresponding ages in mES (left) and HeLa (right) cells. Each dot represents a L1 subfamily and the size of the dot is proportional to its copy numbers around replication origins. Note that some of the youngest L1s with high H3K9me3 eSPAN bias were highlighted in orange. c. Distribution of the length of H3K9me3-bound HO L1s with the lowest H3K9me3 eSPAN bias (red, bottom 25%) and highest H3K9me3 eSPAN bias (blue, top 25%) in mES (left) and HeLa (right) cells. The Y axis was fragmented to better show the details of LINE distribution. Note that L1s with the lowest H3K9me3 bias were shorter than L1s with the highest bias. d. Box plots of H3K9me3 eSPAN bias at HO L1s that were separated into three groups based on their size in mES and HeLa cells. HO L1s were ranked from short to long according to their lengths. The shortest 1/3 was grouped as short (HeLa, n = 3,747; mES, n = 4,225), the middle 1/3 as mid (HeLa, n = 3,767; mES, n = 4,230) and the longest 1/3 as long (HeLa, n = 3,754; mES, n = 4,224). e. Heatmaps of H3K9me3, MPP8 and TASOR ChIP-seq density and H3K9me3 eSPAN bias at HO L1s sorted by L1 length in mES cells, with the size range for long, mid and short L1 groups indicated. The relative position of a full-length L1 was shown in blue. Box plots (d) show the median, 25% and 75% quartiles and minimal and maximal values with p values by two-sided Mann–Whitney–Wilcoxon tests, and Bonferroni correction for multiple comparisons. Each panel is a representative of at least two independent experiments. See Materials and Methods for more details.

Extended Data Figure 5.. Effects of mutating H3K9 methyltransferases on H3K9me3 density and eSPAN bias at L1s in mES cells.

a. SETDB1 depletion reduced H3K9me3 levels dramatically, as detected by immunoblotting. n = 3. b. Heatmaps (left) and average density (right) of H3K9me3 CUT&Tag signals at all the HO L1s in control (shCtr) and SETDB1 knockdown (shSETDB1) mES cells. Heatmaps were sorted by the average H3K9me3 signals of each row in the control sample. The relative position of a full-length L1 was shown in blue. c. Immunoblots of G9a and SUV39h1 to confirm the knockout of G9a, GLP and SUV39h1. Note that antibodies against GLP were not working, but GLP knockout also dramatically reduced the levels of its binding partner, G9a. n = 3. d. Heatmaps (left, sorted as in b) and average density (right) of H3K9me3 CUT&Tag signals at all HO L1s in WT, G9a KO, GLP KO and SUV39h1 KO mES cells. e. Box plots of H3K9me3 density at all HO L1s (n = 12,679) in WT and mutant mES cells. The dashed line indicates the median of H3K9me3 levels in WT cells. f, g. Relative expression of representative repetitive elements in SETDB1 KD (f) and G9a, GLP or SUV39h1 KO (g) mES cells compared to control (shCtr) or WT cells by RT-qPCR analysis. Expression was normalized against shCtr or WT. Data were plotted as mean ± SEM. n = 3-6. h. Immunoblots of ORF1p, the translational products of full-length L1s, in mES cells treated with control or two SETDB1 shRNAs. n = 3. i. H3K9me3 eSPAN bias around replication origins and at HO L1s (n = 12,679) in WT, G9a KO, GLP KO and SUV39h1 KO mES cells. Box plots (e, i) show the median, 25% and 75% quartiles and minimal and maximal values with p values by two-sided Mann–Whitney–Wilcoxon tests, and Bonferroni correction for multiple comparisons. f, g, Two-sided Student’s t test. ****, p < 0.0001. ***, p < 0.001. **, p < 0.01. *, p < 0.05. Each panel is a representative of at least two independent experiments. See Materials and Methods for more details. For gel source data, see Supplementary Figure 1.

Extended Data Figure 6.. TASOR or MPP8 deletion reduced H3K9me3 eSPAN bias in mES cells, while having little effects on H3K9me3 levels.

a. Immunoblots to confirm the knockout of MPP8 and TASOR in mES cells. Note that ORF1p were markedly up-regulated, while total H3K9me3 levels didn’t change to a detectable degree. n = 3. b, c. H3K9me3 CUT&Tag signals at HO L1s in MPP8 KO (b) and TASOR KO (c) mES cells, compared to WT cells. Two repeats for each mutant were shown in the heatmaps (top), with the average H3K9me3 density shown at the bottom. Heatmaps were sorted by the average H3K9me3 signals of each row in WT cells. The relative position of a full-length L1 was shown in blue. d, e. H3K9me3 CUT&Tag signals at HO L1s in MPP8 KO (b) and TASOR KO (c) mES cells, compared to WT cells. Two repeats for each mutant were shown in the heatmaps (top, sorted as in b), with average density shown at the bottom. L1s with reduced H3K9me3 levels for more than 1.5-fold were grouped as Down and those without significant changes were grouped as no-difference (No-diff). Note that less than 130 (~1%) L1s showed significant reduction of H3K9me3 density. f. Snapshots of H3K9me3 CUT&Tag signals and eSPAN bias in WT, MPP8 KO and TASOR KO mES cells at three loci. Note that H3K9me3 CUT&Tag in MPP8 KO and TASOR KO were performed in separate batches with their corresponding H3K9me3 CUT&Tag in WT cells (WT1 and WT2 shown for more accurate comparisons. g. Box plots of H3K9me3 eSPAN bias at two groups of HO L1s (n = 12,679) in WT, MPP8 KO (left) and TASOR KO (right) mES cells. An average of two independent repeats were shown and L1s were grouped as in d, e. Box plots (g) show the median, 25% and 75% quartiles and minimal and maximal values with p values by two-sided Mann–Whitney–Wilcoxon tests, and Bonferroni correction for multiple comparisons. Each panel is a representative of at least two independent experiments. See Materials and Methods for more details. For gel source data, see Supplementary Figure 1.

Extended Data Figure 7.. TASOR or MPP8 depletion reduced H3K9me3 eSPAN bias at HO L1s in HeLa cells.

a. H3K9me3 density at HO L1s in WT and TASOR KO HeLa cells based on published CUT&RUN datasets. Heatmaps (left) were sorted by the average H3K9me3 signals of each row in WT cells, with average density shown at the right. L1s were separated based on the effects TASOR KO on H3K9me3 levels. Of the 393 TASOR regulated H3K9me3 loci identified by Douse et al. using a cutoff of log2 fold-change < −1, we found that 119 HO L1s identified in this study were located at these loci and defined them as the Down group. All other HO L1s were grouped as No-diff. b. H3K9me3 density at HO L1s in WT, MPP8 KO and TASOR KO HeLa cells based on published ChIP-seq datasets, with heatmaps (top, sorted as in a) and average density (bottom) shown. HO L1s were grouped as in a. c. Immunoblots to confirm the knockdown (KD) of MPP8 and TASOR in HeLa cells. Note that while sgRNAs targeting MPP8 or TASOR were used to generate these cells, cells were pooled after selection, instead of cloned. Therefore, MPP8 and TASOR were only depleted and labeled as KD in HeLa cells. Note that ORF1p were markedly up-regulated, while total H3K9me3 levels didn’t change to a detectable degree. n = 3. d. H3K9me3 CUT&Tag signals at HO L1s in WT, MPP8 KD and TASOR KD HeLa cells. The datasets were generated in this study and HO L1s were grouped as in a. e. H3K9me3 CUT&Tag signals at HO L1s separated based on L1 length in WT, MPP8 KD and TASOR KD HeLa cells. HO L1s longer and shorter than the medium length were grouped as long (n = 5,615) and short (n = 5,653), respectively. f. Average H3K9me3 eSPAN bias around all 2,809 replication origins in WT, MPP8 KD and TASOR KD HeLa cells. g, h. Box plots of H3K9me3 eSPAN bias at HO L1s in WT, MPP8 KD and TASOR KD HeLa cells. HO L1s (n = 11,268) were grouped based on the effects of MPP8/TASOR KD on H3K9me3 density defined in a (g), or based on L1 length, as defined in e (h). Box plots (g, h) show the median, 25% and 75% quartiles and minimal and maximal values with p values by two-sided Mann–Whitney–Wilcoxon tests, with Bonferroni correction for multiple comparisons. Each panel is a representative of at least two independent experiments. See Materials and Methods for more details. For gel source data, see Supplementary Figure 1.

Extended Data Figure 8.. The HUSH complex is enriched at the leading strands of DNA replication forks.

a. Detection of the HUSH complex subunits at replication forks based on published iPOND (isolation of proteins on nascent DNA) and NCC (nascent chromatin capture) datasets. Numbers of peptides identified were shown. N.D., not detected. b. Heatmaps of normalized eSPAN density of H3K9me3, MPP8, TASOR and Flag-TASOR at HO L1s, sorted by L1 length. The relative position of a full-length L1 was shown in blue. c. Average MPP8 eSPAN bias around all 1,928 replication origins in mES cells. d-g. Correlations between the biases of TASOR eSPAN and H3K9me3 eSPAN (d, e) or the biases between MPP8 eSPAN and H3K9me3 eSPAN (f, g) in mES cells. Each dot represents a 1 kb bin (d, f) or a HO L1 (e, g) within the 1,928 initiation zones (−100 kb, 100 kb). Spearman’s rank correlation coefficient was shown. p < 2.2e-16. h. Average MPP8 and TASOR eSPAN bias around all 2,809 replication origins in HeLa cells. i. A snapshot of H3K9me3 ChIP-seq and calculated eSPAN biases of H3K9me3, MPP8 and TASOR in HeLa cells. OK-seq bias was shown to mark origin location. j-m. Correlations between the biases of MPP8 eSPAN and H3K9me3 eSPAN (j, k) or between TASOR and H3K9me3 (l, m) in HeLa cells. Each dot represents a 1 kb bin (j, l) or a HO L1(k, m) within the 2,809 initiation zones in HeLa cells (−100 kb, 100 kb). Spearman’s rank correlation coefficient was shown. p < 2.2e-16.

Extended Data Figure 9.. Effects of POLE3 or POLE4 deletion/depletion on H3K9me3 density and H3K9me3 eSPAN bias in mES and HeLa cells.

a. Immunoblots of POLE3 and POLE4 to confirm their deletion in mES and depletion in HeLa cells. Note that cloned ES cells (KO) and pooled HeLa cells (KD) were used for analysis and that H3K9me3 levels remained largely unaffected in the mutant cells. n = 3. b, c. H3K9me3 CUT&Tag (b) or CUT&RUN (c) signals at HO L1s in WT, POLE3 KO and POLE4 KO mES cells. Heatmaps (left) were sorted by the average H3K9me3 signals of each row in WT cells, with average density shown at the bottom. Note that very little changes of H3K9me3 levels were observed in the mutants. d. H3K9me3 CUT&Tag signals at HO L1s in WT, POLE3 KD and POLE4 KD HeLa cells, with the heatmaps (top, sorted as in b) and average density at HO L1s (bottom) shown. HO L1s were grouped as long and short, as defined in Extended Data Figure 7e. e, f. H3K9me3 CUT&Tag signals at HO L1s in WT, POLE3 KO (e) and POLE4 KO (f) mES cells. Heatmaps (left, sorted as in b) and average density (right) were shown. HO L1s were separated into two groups based on the effects POLE3 or POLE4 KO on H3K9me3 levels at HO L1s, with a reduction of more than 1.5-fold defined as the Down group and the rest of L1s within this cutoff being grouped as No-diff group. Note that less than 50 (~0.4%) HO L1s showed a marked reduction of H3K9me3 density and therefore the eSPAN bias was not calculated at this group separately. g. H3K9me3 eSPAN bias around replication origins (top) and at HO L1s (bottom, n = 11,268) in WT, POLE3 KD and POLE4 KD HeLa cells. Long and short HO L1 elements were defined as in Extended Data Figure 7. Box plots show (g) the median, 25% and 75% quartiles and minimal and maximal values with p values by two-sided Mann–Whitney–Wilcoxon tests, and Bonferroni correction for multiple comparisons. Each panel is a representative of at least two independent experiments. See Materials and Methods for more details. For gel source data, see Supplementary Figure 1.

Extended Data Figure 10.. Pol ε coordinates with the HUSH complex for asymmetric H3K9me3 distribution.

a-d. A correlation of the reduction of H3K9me3 eSPAN bias between MPP8 KO (a, b) or TASOR KO (c, d) and POLE3 KO, with each mutant compared to WT mES cells. Each dot represents a 1 kb bin (a, c) or a HO L1(b, d) within the 1,928 initiation zones. Spearman’s rank correlation coefficient was shown. p < 2.2e-16. e-h. Correlation of the reduction of H3K9me3 eSPAN bias between MPP8 KO (e, f) or TASOR KO (g, h) and POLE4 KO compared to WT mES cells. Each dot represents a 1 kb bin (e, g) or a HO L1(f, h) within the 1,928 initiation zones. Spearman’s rank correlation coefficient was shown. p < 2.2e-16. i. H3K9me3 eSPAN bias around replication origins (top) and at HO L1s (bottom, n = 12,679) in WT, POLE3 KO, MPP8 KO and MPP8 KO/POLE3 KO double mutant mES cells. j. Alignment of the protein sequences surrounding an unstructured region of TASOR in 8 different species. M.m., Mus musculus, R.n., Rattus norvegicus, H.s., Homo sapiens, P.t., Pan troglodytes, C.f., Canis familiaris, B.t., Bos taurus, G.g., Gallus gallus, X.l., Xenopus laevis. A predicted alpha helix was indicated and conservancy scores were shown at the bottom. Note that a reported domain that is responsible for the binding of Periphilin-, another HUSH subunit, is in the region. k. Amino acid sequences of the TASOR mutations generated to analyze their effects of mutations on the TASOR-Pol ε binding. The mutated or deleted amino acids were highlighted in red. l, m. TASOR M1/M2 (l) and M3 (m) mutations compromised TASOR interaction with Pol ε subunits, but not MPP8 or PPHLN1. TASOR KO or IgG was used as a negative control. * indicates bands from IgG light or heavy chains. Note that M3 mutation, which contains 18-amino acid deletion in TASOR, caused a major shift of the TASOR band on the gel. n = 3. n. cDNA products of a TASOR fragment amplified from WT or TASOR mutant mES cells for expression TASOR fragments used in the GST pull down assays in Figure 4d. o. Average H3K9me3 eSPAN bias around all 1,928 replication origins in WT or POLE4 KO mES cells treated with triptolide (0.5 μM for 45 min). DMSO was added as a control. p. Average H3K27m3 eSPAN bias around all 1,928 replication origins in MCM2–2A mutant mES cells treated with triptolide (0.5 μM for 45 min). DMSO was added as a control. Note that H3K27me3 eSPAN bias towards the leading strand in MCM2-2A cells is much bigger than that in WT cells, due to the defective transfer of parental histones to the lagging strand, as previously reported^,. q. H3K9me3 eSPAN bias at HO L1s (n = 12,679) in WT or POLE4 KO mES cells treated with triptolide. Note that while triptolide didn’t affect overall H3K9me3 bias around origins, H3K9me3 bias at HO L1s at replicating origins was reduced, suggesting that triptolide affects asymmetric H3K9me3 distribution at selective genomic loci. Box plots (i, q) show the median, 25% and 75% quartiles and minimal and maximal values with p values by two-sided Mann–Whitney–Wilcoxon tests, and Bonferroni correction for multiple comparisons. Each panel is a representative of at least two independent experiments. See Materials and Methods for more details. For gel source data, see Supplementary Figure 1.

Extended Data Figure 11.. Linking asymmetric H3K9me3 segregation at HO Ls to their silencing during S phase.

a. Isolation of mES cells at G1, S and G2 phases based on the expression of Cdt1-mKO2 and Geminin-mAG1 by flow cytometry. Cells were sorted based on the expression of these two cell cycle indicators. G1 phase cells only express Cdt1, but not Geminin. S phase cells express medium levels of Geminin, but not Cdt1, and G2 phase cells express the highest levels of Geminin. To increase the purity of S and G2 phase cells, we used a stringent gating strategy as shown with isolated G1, S and G2 phase of cells accounting for ~10%, ~15% and ~10% of total cells, respectively. b. Heatmaps of differentially expressed L1 elements in MPP8 KO, POLE3 KO and POLE4 KO versus WT mES cells based on GRO-seq analysis. Numbers of total differentially expressed L1s in each mutant were shown. Please note that all L1 elements were used in the analysis. c. Relative expression of representative repetitive elements in MPP8 KO, MPP8 W80A or TASOR M3 mutant cells compared to WT mES cells by RT-qPCR. The expression was normalized against WT. Data were plotted as mean ± SEM. n = 5-9. d. Snapshots of H3K9me3 signals (both ChIP-seq and CUT&Tag), H3K9me3 eSPAN bias and GRO-seq signals at three L1 elements in WT or mutant mES cells. The three up-regulated L1s were highlighted. e. Relative expression of representative repetitive elements in POLE3 KD and POLE4 KD cells compared to WT HeLa cells detected by RT-qPCR. The expression was normalized against WT. Data were plotted as mean ± SEM. n = 3. f. The expression of HO L1s (n = 2,662) in WT, POLE3 KO and POLE4 KO mES cells detected by GRO-seq after excluding the ones located within the transcribed regions of up-regulated genes in POLE3 KO or POLE4 KO mES cells defined by RNA-seq. g. The expression of HO L1s (n = 2,309) in WT, POLE3 KO and POLE4 KO mES cells detected by GRO-seq after excluding the ones located within any actively transcribed genes from analysis (cutoff: TPM > 0.5). h. The expression of full-length (≥ 6 kb) HO L1s (n = 703) with their own promoters in WT, POLE3 KO and POLE4 KO mES cells. See Materials and Methods section for more details for panels f-h. c, e, Two-sided Student’s t test. ****, p < 0.0001. ***, p < 0.001. **, p < 0.01. *, p < 0.05. Box plots (f-h) show the median, 25% and 75% quartiles and minimal and maximal values with p values by two-sided Mann–Whitney–Wilcoxon tests, and Bonferroni correction for multiple comparisons. Each panel is a representative of at least two independent experiments. See Materials and Methods for more details.

Extended Data Figure 12.. Effects of POLE3 KO, POLE4 KO and TASOR/MPP8 mutants on L1 expression and retrotransposition.

a. The expression of HO L1s (n = 2,681) in WT, MPP8 KO and MPP8 KO/POLE3 KO double mutant mES cells detected by GRO-seq. b. Overlaps between the up-regulated HO L1s in MPP8 KO and POLE3 KO (top) or between MPP8 KO and POLE4 KO (bottom) mES cells detected by GRO-seq. P values by hypergeometric test. c-e. Comparison of properties (L1 length, TASOR density, and H3K9me3 eSPAN bias) of HO L1s whose expression is up-regulated in MPP8 KO (c), POLE3 KO (d) or POLE4 KO (e) mES cells to those HO L1s without changes in expression in the corresponding mutants. L1s with more than 1.5-fold increase in expression were grouped as Up and those within the 1.5-fold threshold were grouped as No-diff. f. Relative expression of HO L1s (n = 2,681) in MPP8 KO versus WT mES cells at G1, S or G2 phase of the cell cycle detected by GRO-seq. The dashed line indicates no changes compared to WT cells (0). g. Snapshots of GRO-seq signals at the indicated L1 elements in G1, S and G2 phases of WT, POLE4 KO and MPP8 KO mES cells. H3K9me3 ChIP-seq signals and eSPAN biases at these two loci in both WT and mutant cells were also shown. h. Relative L1 mobility in WT, POLE3 KD and POLE4 KD HeLa cells as measured by dual-luciferase reporter assays. Data were plotted as mean ± SEM. n = 8. i. The H3K9me3 eSPAN bias correlates with L1 integration at the leading strands. Absolute values of H3K9me3 eSPAN bias were separated into eleven equal intervals from 0 to 1 (X axis). The fraction of insertions where (+) strand of L1 cDNA integrated into the predominant leading strand template (Y axis) was plotted at each of the matching H3K9me3 bias interval. j. Overlaid violin plots of H3K9me3 eSPAN bias frequency distributions for L1 integrations into the reference genome. Observed L1 insertions in HeLa cells were stratified by the integration strand. The colored lines identify L1 integration into the Watson (orange) and Crick (green) strands of human genome, which means that L1 endonuclease cleaved the opposite strands, i.e., the Crick and Watson strands, respectively. All violin plots were adjusted to have the same total area and vertical lines denote the distribution medians. k. Relative γ-H2AX signal intensity measured by immunofluorescence in WT (n = 551), POLE3 KD (n = 286) and POLE4 KD (n = 362) HeLa cells. WT cells were treated with 1 mM hydroxyurea (HU, n = 361) for 1 h as a positive control. Data were plotted as mean ± SD. Box plots show the median, 25% and 75% quartiles and minimal and maximal values. a, c-f, k, p values by two-sided Mann–Whitney–Wilcoxon tests, and Bonferroni correction for multiple comparisons. h, Two-sided Student’s t test. ****, p < 0.0001. **, p < 0.01. Each panel is a representative of at least two independent experiments. See Materials and Methods for more details.

Figure 1.. H3K9me3 is transferred preferentially to the leading strands.

a. Left: a schematic diagram of the eSPAN procedure and calculation of eSPAN bias. In this hypothetical model, parental H3K9me3 is transferred to leading strands of DNA replication forks, with two nucleosomes at each leading or lagging strand drawn for simplicity. W and C: sequence reads of Watson and Crick strands, respectively. Right: average bias of H3K9me3, H3K9me2, H3K27me3 and H4K20me3 eSPAN surrounding 1,928 replication origins (−100 kb to 100 kb), with two independent repeats (blue and red) shown. b. Heatmaps of eSPAN bias for H3K9me3, H3K9me2, H3K27me3 and H4K20me3 centered around each of the 1,928 replication origins sorted based on replication efficiency defined by OK-seq (right). c. A snapshot of H3K9me3 and H3K27me3 ChIP-seq signals and calculated H3K9me3 and H3K27me3 eSPAN bias, with OK-seq bias used to indicate origin location and DNA replication direction (shown by arrow). L1 elements (≥ 1 kb) at this locus were shown at the bottom with their transcription direction indicated. d. Average bias of H3K9me3, H3K27me3 and H4K20me2 eSPAN signals in both HeLa cells and primary mouse B cells surrounding 2,809 and 1,073 replication origins (−100 kb to 100 kb), respectively.

Figure 2.. H3K9me3 asymmetry occurs in S phase and at LINE elements.

a. Average H3K9me3 eSPAN bias around origins in mES cells at different cell cycle stages. The respective cell cycle profiles were shown on the bottom. b. Enrichment of different repetitive elements at 1 kb bins with high (top quartile, 25%) and low (bottom quartile, 25%) H3K9me3 eSPAN bias in mES cells. Fold enrichment is defined as the ratio between calculated and expected enrichment. com., complexity; rep., repeats; LTR, long-terminal repeat; SINE, short interspersed nuclear element; DNA, DNA transposon; scRNA, small cytoplasmic RNA; snRNA, small nuclear RNA; tRNA, transfer RNA; rRNA, ribosomal RNA; srpRNA, signal recognition particle RNA; RC, rolling circle; unknown, repeats without a recognizable TE signature. c. Box plots of H3K9me3 eSPAN bias at replicated LINEs (n = 20,513) and all other TEs (n = 71,824) that were separated into co-direction (CO) and head-on (HO) groups based on the transcription direction of each TE unit and replication fork direction. d. Box plots of L1 length, TASOR and H3K9me3 ChIP-seq density, and TASOR eSPAN bias for the top and bottom quartile of HO L1s with the highest and lowest eSPAN bias (total HO L1s = 12,679). e. Box plots of H3K9me3 eSPAN bias at HO L1s with high TASOR density (Q4, n = 3,170) and low TASOR density (Q1, n = 3,170) based on TASOR ChIP-seq. f. Heatmaps of H3K9me3, MPP8 and TASOR ChIP-seq density, and H3K9me3 eSPAN bias at HO L1s sorted by TASOR levels. The relative position of a full-length L1 was shown in blue. Box plots (c-e) show the median, 25% and 75% quartiles and minimal and maximal values with p values by two-sided Mann–Whitney–Wilcoxon tests. c, Bonferroni correction for multiple comparisons. Each panel is a representative of at least two independent experiments. See Materials and Methods for more details.

Figure 3.. The HUSH complex regulates H3K9me3 asymmetry.

a. H3K9me3 eSPAN bias around replication origins (left) and at HO L1s (right, n = 12,679) in control (shCtr) or SETDB1 knockdown (shSETDB1) mES cells. b. H3K9me3 eSPAN bias around replication origins (left) and at HO L1s (right, n = 12,679) in WT, TASOR and MPP8 mutant mES cells. c. Heatmaps of normalized eSPAN density of H3K9me3, MPP8, TASOR and Flag-TASOR at HO L1s, sorted by TASOR eSPAN density in mES cells. The relative position of a full-length L1 was shown in blue. d. TASOR eSPAN bias around replication origins (left) and at HO L1s (right, n = 12,679) with high and low TASOR ChIP-seq density defined as in Figure 2e, in WT mES cells. e. A snapshot of H3K9me3, TASOR and MPP8 ChIP-seq signals and calculated eSPAN bias at the indicated locus, with OK-seq bias marking origin location. All L1s at this locus were shown at the bottom. Box plots (a, b, d) show the median, 25% and 75% quartiles and minimal and maximal values with p values by two-sided Mann–Whitney–Wilcoxon tests. b, Bonferroni correction for multiple comparisons. Each panel is a representative of at least two independent experiments. See Materials and Methods for more details.

Figure 4.. Pol ε coordinates with the HUSH complex for asymmetric H3K9me3 transfer.

a. H3K9me3 eSPAN bias around replication origins (left) and at HO L1s (right, n = 12,679) in WT, POLE3 KO and POLE4 KO mES cells. b. Interactions between MPP8 and Pol ε determined by co-immunoprecipitations. Anti-Rabbit IgG was used as a negative control. * indicates the IgG light chain. n = 3. c. TASOR KO or MPP8 chromodomain mutation (W80A) compromised the MPP8-Pol ε interaction. MPP8 KO was used as a negative control. * indicates the IgG light chain. n = 3. d. TASOR site-specific mutations (M1, M2, and M3) compromised the TASOR-Pol ε interaction as determined by in vitro GST pull-down assays. GST-Reg α was used as a negative control. n = 5. e, f. H3K9me3 eSPAN bias around replication origins (left) and at HO L1s (right, n = 12,679) in WT and MPP8 W80A (e) or three TASOR mutants (f) mES cells. Box plots (a, e, f) show the median, 25% and 75% quartiles and minimal and maximal values with p values by two-sided Mann–Whitney–Wilcoxon tests. a, f, Bonferroni correction for multiple comparisons. Each panel is a representative of at least two independent experiments. See Materials and Methods for more details. For gel source data, see Supplementary Figure 1.

Figure 5.. Effects of asymmetric H3K9me3 distribution on L1 silencing and retrotransposition.

a. Expression of HO L1s (n = 2,681) in WT mES cells at G1, S or G2 phase of the cell cycle detected by GRO-seq. Note that HO L1s with expression lower than the cutoff (CPM = 0.1) were excluded from analysis. The same number of L1s was used for analysis in panels b, d. b. Expression of HO L1s (n = 2,681) in WT and mutant mES cells detected by GRO-seq. c. Relative expression of representative repetitive elements of different classes in WT and mutant mES cells by RT-qPCR. The expression was normalized against WT (mean ± SEM. n = 3). Sat., satellite. mLINE1, mouse L1 elements. d. Relative expression of HO L1s (n = 2,681) in POLE4 KO versus WT mES cells at G1, S or G2 phase of the cell cycle. The dashed line indicates no changes compared to WT (0). e. H3K9me3 eSPAN bias changes were negatively correlated with L1 activation in the mutant cells. Each dot represents an expressed HO L1 within the 1,928 initiation zones. f. Relative L1 mobility in WT and mutant mES cells (mean ± SEM. n = 6). g. Schematic working model. In WT cells, asymmetric H3K9me3 distribution to head-on L1s at the leading strand is regulated by low level transcription, the HUSH complex and Pol ε, which is important to inhibit L1 transcription and retrotransposition during S phase. Box plots (a, b, d) show the median, 25% and 75% quartiles and minimal and maximal values with p values by two-sided Mann–Whitney–Wilcoxon tests, and Bonferroni correction for multiple comparisons. Each panel is a representative of at least two independent experiments. c, f, Two-sided Student’s t test. ****, p < 0.0001. ***, p < 0.001. **, p < 0.01. *, p < 0.05. n.s., not significant, p > 0.05. See Materials and Methods for more details.