LINE-1 transcription activates long-range gene expression

a, Heatmap showing the relative expression of individual full-length L1s across three indicated human cell lines. Z-score transformed transcripts per million (TPM) values were used to represent L1 expression. b, L1 RNA quantification in K562 using the strand-specific polyA+ RNA-seq. Data are representative of two independent experiments. The center value is the mean. P, two-tailed t test. n.s., non-significant. c, Fluorescence-activated cell sorting results show that over 99% of the L1-GFP K562 cell clones express blue fluorescent protein (BFP)-Cas9 and L1-GFP. d, Chromosome maps showing the genomic distribution of 51 L1-GFP reporter insertions in the K562 screening cell line, measured by long-read sequencing. The individual genomic locations are listed in Supplementary Table 1. e, Normalized L1-GFP expression (left) and line plots showing GFP intensity (right) in control and mutant K562 as indicated. The center value is the mean. n = 2 biological replicates for MPP8, n = 3 biological replicates for SAFB and control. Knockout of known L1 suppressors (SAFB, MPP8) increased L1 expression. f, Correlation of the sgRNA counts of each sample in K562 genome-wide screens, relative to Fig. 1a. g, Enrichment values of the top 80 L1 suppressors (blue) and the top 80 activators (red) were calculated from two independent K562 genome-wide screens with 4 independent sgRNAs per gene. Diamond represents the mean value; whisker represents the fold range for sgRNAs. Previously characterized L1 regulators are marked with stars. h, Cytoscape for visualizing protein–protein interactions for L1 suppressors (blue) and L1 activators (red). i,j, Venn diagram showing the number of candidate L1 activators (i) and suppressors (j) identified in this screen against L1 expression (blue), and that identified from one previous screen against L1 retrotransposition¹⁶ (yellow) in K562. k,l, Gene ontology (biological process) enrichment analysis of the unique activators (k) and unique suppressors (l) identified in this screen against L1 expression, and that identified from one previous screen against L1 retrotransposition in K562 (ref. ¹⁶). P, one-tailed hypergeometric test.

a, Bar plots showing the percentage of unique-mapped and multiple-mapped reads in indicated full-length L1PA subfamilies. b, Box plot showing the percentage of unique-mapped reads in individual full-length L1PAs (left), and the numbers of all elements and those that surpass the 5-read-count cutoff within each L1 subfamily (right). Center lines represent the median value, box limits represent the 25th and 75th percentiles and whiskers denote minima and maxima (1.5× the interquartile range). n = 2 biological replicates. c, Representative genome browser snapshots of the RNA-seq reads alignment over L1s. The color scale indicates the mapping quality score MAPQ for each read pair. MAPQ = −10 log₁₀(p), where p is the probability that true alignment belongs elsewhere. d, Heatmaps showing expression change of significantly changed L1PAs using unique mapping (top) or multiple mapping (bottom) in indicated mutant K562. n = 2 biological replicates per gene. P < 0.05, DESeq2 analysis. e, Mean splicing efficiency of the RNA-seq samples, calculated using SPLICE-q⁷⁹. f, Intronic read fraction of the RNA-seq samples, calculated using RNA-SeQC⁸⁰. Together with e, these results indicated the high quality of our RNA-seq library. g, Bubble plots showing expression change of L1, LTR and SINE subfamilies in indicated mutant K562. Multi-mapped reads were used to quantify transposons at the subfamily level. Colored circles represent the log₂ fold change, with areas proportional to P value (DESeq2 analysis). LTR, long terminal repeat; SINE, short interspersed nuclear element. n = 2 or 3 biological replicates per gene. h, Heatmap showing expression change at individual L1 locus in indicated mutant K562 (labeled and aligned same as in d). Hierarchical clustering was performed by ‘ward.D2’ method. RNA-seq reads were unique-mapped to individual L1 locus. Only significantly changed full-length L1s in at least one KO experiment were included (n = 222). P < 0.05, DESeq2 analysis. i,j, RNA-seq genome browser snapshots at six L1 loci from control and indicated mutant K562, illustrating that candidate L1 suppressors (i) and activators (j) selectively and cooperatively control L1 expression. The experiment was repeated once with similar results.

a, Western blotting showing the level of DBR1 and endogenous L1 ORF1p in the K562 cells infected with two independent Cas9/sgRNAs targeting DBR1 and two independent non-target control Cas9/sgRNAs, independently repeated once with similar results. HSP90 is loaded as control. b, Donut chart showing that DBR1-regulated L1s (P < 0.05, DESeq2 analysis) are predominantly localized in introns in K562. c, Size distribution of DBR1-regulated L1s (P < 0.05, DESeq2 analysis) and all genomic L1s (based on hg38 RepeatMasker annotation) in K562. P value, two-tailed Kolmogorov–Smirnov test. d, The read counts of total lariats and L1-containing lariats, obtained from the ribosome-depleted RNA-seq data by a lariat read-count analysis strategy⁸¹ and normalized against unmapped singlets in control and DBR1 knockout K562 cells. n = 2 biological replicates. The center value is the mean. e, A working model that DBR1 restricts intronic L1 expression by debranching lariats into linear RNAs that are subsequently degraded. In DBR1 mutant cells, the level of L1-containing lariats is elevated, and the L1 within the lariat is proposed to be translated into L1 proteins. f, Donut chart showing genomic distribution of the DIS3-restricted L1s (P < 0.05, DESeq2 analysis) in K562 cells. g, Density plot showing size distribution of the DIS3-regulated L1s (P value < 0.05, DESeq2 analysis) and all genomic L1s (based on hg38 repeatmasker annotation) in K562 cells. P value, two-tailed Kolmogorov–Smirnov test. h, Schematic of DIS3 domains and the two point mutations to inactivate the PIN and RNB domains, respectively²¹. i, RNA-seq showing the relative expression level of full-length L1s in DIS3 knockdown HEK293 cells, rescued with the intact DIS3 or DIS3 mutants with inactive PIN or RNB domains, as indicated. RNA-seq reads were unique-mapped to individual L1 locus. As indicated, the PIN domain with endonucleolytic activity accounts for the L1 restriction. j, Aggregated line plot showing average enrichment of the DIS3 PAR-CLIP signal over full-length L1s. Genome-wide mappability displays genome-wide 100-bp-short-read alignability (data from UCSC). The results indicate that DIS3 directly binds L1 RNAs.

Source data

a,b, Dot plot/box plot of STARR-seq read counts (left) and enrichment (right) at the L1 subfamily (a) and L1 individual (b) level using K562 STARR-seq datasets (2 PE100 and 1 PE150). The STARR-seq reads were unique-mapped to individual L1 locus. Center lines represent the median value, box limits represent the 25th and 75th percentiles and whiskers denote minima and maxima (1.5× the interquartile range). n = 2 biological replicates. CPM, read counts per million mapped reads. c, Scatter plot showing counts of STARR-seq peaks and −log₁₀(P value) for all repeats at the subfamily level in indicated cell lines. Dots, repeat subfamilies; large dots, L1 subfamilies; color, repeat classes. P, hypergeometric test. d, Aggregated line plots showing STARR-seq signals over full-length L1s at the subfamily level in indicated cell lines. e, Heatmaps of STARR-seq signals over full-length L1s in indicated cell lines, sorted and aligned as in Fig. 2d. f, Heatmap showing the proportion of full-length and non-full-length L1s that exhibit STARR-seq signals in six indicated human cell lines. The L1 subfamilies are classified into two groups according to their length: longer than 6 kb (full-length) or shorter (non-full-length). While nearly all full-length L1s contain L1 5′ UTR and overlap with STARR-seq signals, some non-full-length L1s (such as L1P2 and L1PBa1) also contain L1 5′ UTR and overlap with STARR-seq peaks. g, Enrichment of DNA binding motifs at the full-length L1s with STARR-seq signals in K562, with the full-length L1s not overlapped with STARR-seq peaks as control. Plotted is the log₂ transformed fold enrichment and log₁₀ transformed P value of each protein binding motif (hypergeometric test). The enriched motifs for transcription factors are specified. Zf, zinc finger. h, Analysis of the ChIP-Atlas data showing that transcription factors bind full-length L1s. Colored circles represent log₂ fold enrichment, with areas proportional to Q value (hypergeometric test). For a single protein, log₂ fold enrichment is calculated using the median value from different cells or experiments. Q-values are adjusted with the Benjamini–Hochberg method.

a, Western blotting showing that CTBP1 KO increases L1 protein expression in K562. b, FACS showing CTBP1 KO increases L1-GFP expression in K562. Center line, mean. n = 2 biological replicates. c, qRT–PCR data showing normalized expression of L1Hs, L1PA15-16, L1PB and L1M subfamilies in control (n = 3) and CTBP1 KO cells (n = 3) with or without reverse transcription (RT). P, two-tailed t test. d, Box plot showing that CTBP1 KO mainly derepresses L1PA1–3 among all repeats in K562. Box plots show median and interquartile range (IQR), and whiskers are 1.5-fold IQR. e, Density plot showing size distribution of the CTBP1-regulated L1s (P < 0.05, DESeq2 analysis) and other L1s in K562. P, two-tailed Kolmogorov–Smirnov test. f, Aggregated line plots showing indicated ChIP enrichments over the CTBP1-regulated L1PAs (P < 0.05, DESeq2 analysis) and all full-length (FL) L1s in K562. g, CTBP1 KO increases expression of genes surrounding the derepressed L1s. Plots show median and interquartile range (IQR) of gene expression change alongside derepressed L1s (black box plots), compared with background (red). h, Bar plot showing the number of full-length L1s detected by H3K27ac ChIP–seq and polyA+ RNA-seq (TPM > 0.05) in K562. i, Zoom-in views of the RNA-seq genome browser tracks at NRN1, its cognate distal L1PA2 and one nearby FARS2 locus, in control and CTBP1 KO K562, independently repeated once with similar results. j, RT–qPCR showing normalized expression of NRN1 and its nearby FARS2 in control (n = 3) and CTBP1 KO cells (n = 3). P, two-tailed t test. n.s., not significant. k,l, Deleting the L1PA2 element, as shown in i and confirmed by the gel image (k), decreases NRN1 expression in NCCIT (l). n = 4 biological replicates × 3 technical replicates. ***P < 0.001; unpaired two-tailed t-test. Error bars, standard deviation (s.d.). m,n, Representative Hi-C interaction matrix and ChIP–seq IGV tracks around NRN1-L1 in wild-type (WT) condition (m) and Sanger sequencing results showing the L1 deletion sites with sgRNA target sequences in red (n).

Source data

a, Western blotting showing KMT2D KO increases L1 protein expression in K562. b, FACS showing KMT2D KO increases L1-GFP expression in K562. Center line, mean. n = 2 biological replicates. c, Density plot showing size distribution of the KMT2D-regulated L1s (P < 0.05, DESeq2 analysis) and other L1s in K562. P, two-tailed Kolmogorov–Smirnov test. d, Bubble plot showing expression changes of different L1 subfamilies in mouse MA1L cells. Data from ref. ⁸². e, Heatmaps showing indicated ChIP–seq enrichments in control and KMT2D KO K562 cells, sorted by the increase of H3K4me1 signals upon KMT2D KO and centered on the H3K4me1 peaks. f, Heatmaps showing H3K4me1 ChIP–seq enrichment over transposable elements, centered and sorted as in e. The yellow rectangle indicates the full-length L1s with increased H3K4me1 peaks upon KMT2D knockout in K562. g, Heatmaps showing indicated ChIP enrichments in control and KMT2D KO K562 over full-length L1s, sorted by WT H3K4me1 ChIP and aligned at the L1 5′ end. h, Frequency histogram of absolute distances from each L1 to the nearest gene with increased (red) and decreased expression (blue) upon KMT2D KO in K562. Background expectations of genome-wide L1 distribution to the nearest upregulated (green) and downregulated (black) gene upon KMT2D KO are shown, with dashed line indicating 95% confidence interval. P, two-tailed Kolmogorov–Smirnov test; n.s, non-significant. i, Box plots showing expression change for genes with an L1 which shows increased H3K27ac ChIP signals upon KMT2D KO in K562, within indicated genomic distance. Plots show median and interquartile range (IQR), and whiskers are 1.5-fold IQR. Related statistical significance analysis is in Supplementary Table 5. j, Bubble plots showing expression changes of different L1 subfamilies in mESC. Data from refs. ^83,84,85. k, Scatter plot showing METTL3/14 ChIP–seq enrichment over repeats in MCF-7 cells. Dots represent repeat subfamilies, colored by repeat classes. P, hypergeometric test. l, Box plots showing expression changes for genes with an L1 which shows decreased expression upon METTL3 KO in K562, within indicated genomic distance. Plots show median and interquartile range (IQR), and whiskers are 1.5-fold IQR. Related statistical significance analysis is in Supplementary Table 5.

Source data

a, Schematic of triple sgRNAs targeting L1Hs 5′ UTR, used in CRISPRa/i. b, RT–qPCR results showing L1Hs expression upon CRISPRa/i in NCCIT. P, unpaired two-tailed t-test. Data represent median ± SEM. n = 2 biological replicates × 3 technical replicates. RT, reverse transcription. NC, non-target control. c, Bar plots showing counts of indicated L1s (blue + yellow + gray), counts of L1s targeted by triple sgRNAs (blue + yellow) and counts of L1s with increased H3K27ac peaks upon L1 CRISPRa (blue) in NCCIT. d, Bar plots showing counts of indicated L1s with H3K27ac peaks in NCCIT (blue + yellow + gray), counts of L1s targeted by triple sgRNAs (blue + yellow) and counts of L1s with decreased H3K27ac peaks upon L1 CRISPRi (blue). e, Gene expression change upon L1 CRISPRa/i in NCCIT. Red dots, genes whose expression increased upon CRISPRa (padj < 0.1, DESeq2 analysis) and decreased or remained silenced upon CRISPRi (termed as the L1-regulated genes). f, Plotted as in e, with genes separated into deciles by distance from the nearest L1 with increased H3K27ac peaks upon L1 CRISPRa. g, Box plot showing gene expression change in different deciles, binned as in f. n = 2 biological replicates. Center lines represent the median value, box limits represent the 25th and 75th percentiles and whiskers denote minima and maxima (1.5× the interquartile range). h,i, Representative examples of the L1-initiated chimeric transcripts in NCCIT. j, Box plots showing expression changes for the genes with Mll3/4 binding sites within indicated genomic distances in mESCs. Plots show median and interquartile range (IQR), and whiskers are 1.5-fold IQR. This suggests that the enhancer activity of L1 (shown in Fig. 3d,e) is comparable to that of canonical enhancers regulated by Mll3/4. k, RNA-seq browser tracks at NRN1 and its cognate L1PA2 locus. l, RNA-seq browser tracks at SESN3, its cognate L1Hs and nearby gene locus. m, RT–qPCR results show that CRISPRa targeting the SESN3-L1, rather than one adjacent locus, increases SESN3 expression. Data represent median ± SEM. n = 2 biological replicates × 3 technical replicates. P, two-tailed t-test. ****P < 0.0001. n, IGV tracks in NCCIT after CRISPRa against L1 (data from this study) and LTR5Hs (result from ref. ²⁹). This suggests the importance of L1 transcription, rather than the adjacent LTR5Hs, in the long-range gene activation of HMGCS1.

Source data

a,b, Scatter plot comparing genome-wide compartment scores upon L1 CRISPRa (a) and CRISPRi (b) in NCCIT. NC, non-target control. c, Stacked bar plot showing compartment changes upon L1 CRISPRa/i in NCCIT cells. d, Distribution of the distance between L1s and their looped cognate genes in NCCIT. e, Pie-chart representing number and proportion of L1s with the same or opposite orientation (i), and upstream or downstream (ii) of their looped genes identified in L1 CRISPRa (up) and CRISPRi (down) NCCIT. f, Enrichment analysis of transcription factor binding motifs at the 5′ UTRs of the active L1s associated with genes compared with those that do not. g,h, Box plot showing expression change of L1-looped genes compared to non-looped genes within 2 Mb of the H3K27ac- increased or decreased L1s in L1 CRISPRa (g) and CRISPRi (h), respectively. Plot shows the median and interquartile range (IQR), and whiskers are 1.5 × IQR. P, two-tailed Mann–Whitney–Wilcoxon test. i,j, Scatter plot showing correlation between the fold change of Pol II ChIP–seq signal at individual L1 and the expression fold change of the L1-paired genes in NCCIT upon L1 CRISPRa (i) and CRISPRi (j). k,l, Aggregated line plots showing ChIP enrichments of Pol II (k) and RAD21 (l) over the full-length L1s that contact their cognate gene promoters (left), and the respective overall ChIP–seq peaks (right) in L1 CRISPRa (blue) and non-target control (NC, gray) NCCIT cells. m, Bar plot showing the ChIP–seq enrichment of transcription factors and histone marks at the L1s that make physical contacts with their cognate genes. ChIP–seq datasets in the human embryonic stem cells (H9 hESC) were from ENCODE. n, Hallmark enrichment analysis for the L1-regulated ZGA genes. P, hypergeometric test. o,p, Scatter plot showing correlation between the H3K27ac fold changes at individual L1 locus and the fold changes of paired gene expression in NCCIT upon L1 CRISPRa (o) and CRISPRi (p).

a, Scatter plot showing STARR-seq peak counts and −log₁₀(P) for all repeats in mESCs. Dots, repeat subfamilies; large dots, L1 subfamilies; color, repeat classes. P, hypergeometric test. b, Aggregated contacts between L1s and their looped genes at indicated mouse embryos using the Hi-C data⁴⁹. c, Box plot quantifying the L1–gene contacts as described in b. Plots show median and interquartile range (IQR), and whiskers are 1.5-fold IQR. **P < 0.01, two-tailed t-test. d, Bubble plot showing overlap between the L1-contacted genes (c) and the list of 2C genes⁴³ and DBTMEE v2 transcriptome categories⁸⁶. P, hypergeometric test. e, Volcano plots showing repeat expression changes upon L1 perturbation in mESCs. f, Aggregated Hi-C contacts between L1s and their looped genes and all looping events in control and L1 CRISPRi mESCs. g, Bubble plot showing overlap between all DEGs (differentially expressed genes) by indicated L1 perturbation with the list of 2C genes⁴³ and DBTMEE v2 transcriptome categories⁸⁶. P, hypergeometric test. h, Gene set enrichment analysis (GSEA) among all Hallmark gene sets upon indicated L1 perturbation. NES, normalized enrichment scores. i, Gene set enrichment analysis (GSEA) showing that 2C-specific genes (defined in ref. ⁴³) enrichment in mESCs after indicated L1 perturbation. j, Heatmap showing expression changes of 2C-specific genes and MERVLs upon indicated L1 perturbation. MERVL counted at subfamily level. k, Representative bar plots showing gene expression changes for indicated 2C genes upon indicated L1 perturbation. Data are representative of two or three independent experiments and are median ± SEM. *P < 0.1, **P < 0.01, ***P < 0.001, ****P < 0.0001, one-tailed t test. l, Box plot showing distance distribution of 2C genes sets to the nearest full-length L1 elements. Plots show median and interquartile range (IQR), and whiskers are 1.5-fold IQR. P, two-tailed Mann–Whitney–Wilcoxon test. m, Deleting an L1MdA element, as shown in Fig. 5f and validated by gel image (left), decreases Zscan4a/4c/ps-1 expression. n = 2 biological replicates × 3 technical replicates. ***P < 0.001; unpaired two-tailed t-test. Error bars, standard deviation (s.d.).

a, Gene set enrichment analysis (GSEA) showing that 2C-specific genes were significantly downregulated in the 2-cell mouse embryos after L1 ASO treatment. NES, normalized enrichment score. b, Scatter diagram and box plot (inside) showing the relationship between the distance of gene transcription start site (TSS) to the nearest full-length L1s and the gene expression change upon L1 ASO treatment in mouse 2C embryos. The two vertical lines indicate ±512 kb, respectively. The top density plot shows the distance distribution of the expressed 2C-specific genes (red) and other genes (gray). The box plot shows the median and interquartile range (IQR) of gene expression change, and whiskers are 1.5× IQR. P, two-tailed Mann–Whitney–Wilcoxon test. The 2C-specific genes are significantly enriched within ±512 kb of the full-length L1s, and 96 of them are significantly downregulated upon L1 ASO in mouse 2C embryos. c, Examples of gene clusters of 2C-specific genes, with the derepressed ones in L1-ASO 2-cell mouse embryos labeled in red. d, Chromosome maps illustrating the genomic distribution of the 2C-specific genes, with the derepressed ones in L1-ASO 2-cell mouse embryos labeled in red. e, Representative phase-contrast images of 2.5 and 3.5 dpc mouse embryos in uninjected, control ASO and L1 ASO conditions, from 3 independent experiments. Scale bars, 100 µm. f, Percentage of mouse embryos at 4.5 dpc in SCR (n = 45), RC-L1 ASO (n = 44) and L1 ASO (n = 39). ***P < 0.0001, chi-square test. n.s., not significant. SCR ASO, scramble non-targeting ASO. RC-L1 ASO, the ASO with reverse complementary sequence to the L1 ASO sequence. g, MA plots of expression change of the 2C-specific genes in 2-cell mouse embryos after microinjection of RC-L1 ASO (left) and L1 ASO (right) into zygote, relative to the SCR ASO control. h, Model showing that L1 knockdown causes ZGA defects and development arrest at 2–4 cell stage before blastocyst formation in mouse early embryos. The figure is created with BioRender.com.

Source data