Nuclear RNA catabolism controls endogenous retroviruses, gene expression asymmetry, and dedifferentiation

Abstract

Endogenous retroviruses (ERVs) are remnants of ancient parasitic infections and comprise sizable portions of most genomes. Although epigenetic mechanisms silence most ERVs by generating a repressive environment that prevents their expression (heterochromatin), little is known about mechanisms silencing ERVs residing in open regions of the genome (euchromatin). This is particularly important during embryonic development, where induction and repression of distinct classes of ERVs occur in short temporal windows. Here, we demonstrate that transcription-associated RNA degradation by the nuclear RNA exosome and Integrator is a regulatory mechanism that controls the productive transcription of most genes and many ERVs involved in preimplantation development. Disrupting nuclear RNA catabolism promotes dedifferentiation to a totipotent-like state characterized by defects in RNAPII elongation and decreased expression of long genes (gene-length asymmetry). Our results indicate that RNA catabolism is a core regulatory module of gene networks that safeguards RNAPII activity, ERV expression, cell identity, and developmental potency.

Keywords: 2CLC; Integrator; MERVL; RNA catabolism; elongation; endogenous retrovirus; non-coding RNA; stem cell; totipotent-like cells; transcription-associated RNA degradation.

Figure 1.. Loss of Exosc3 upregulates LTR-containing TEs in mESCs

(A) Volcano plot of differentially expressed genes and TEs between WT and Exosc3 cKO mESCs. Genes that are significantly upregulated and downregulated in Exosc3 cKO mESCs are shown in red and blue, respectively (|log₂-fold change| > 1, p < 0.05, Benjamini-Hochberg correction, calculated using DESeq2). Total number of significantly upregulated and downregulated genes is indicated. (B) Bar plot displaying the number of TEs significantly upregulated in Exosc3 cKO mESCs. (C) Circos plot of PacBio de novo sequenced genome of Exosc3 COIN mESCs compared with mm10. Green represents sequences within the new assembly, but not in mm10 (insertions). Red represents sequences not found in the new assembly but found in mm10 (deletions). The y axis indicates the size of the structural variants in base pairs (positive values for insertions, negative values for deletions). (D) Number of novel inserted (green) and deleted (red) TE copies in the de novo-sequenced genome compared with mm10. (E) Box plots displaying DNA methylation at genes and TEs in WT and Exosc3 cKO mESCs. (F) Characterization of the basal chromatin state of TEs upregulated in Exosc3 cKO using ChIP-seq data. (G) Bar plots displaying RNAPII levels at TEs in WT and Exosc3 cKO mESCs (log₂-fold change of ChIP-seq signal vs. background, 8WG16 antibody, top 15 TEs displayed). p-values were calculated using an unpaired, two-sided Wilcoxon rank-sum test and adjusted using the Benjamini-Hochberg method (**p < 0.01, ****p < 0.0001). (H) Enrichment of RNAPII ChIP-seq at full-length MERVL elements (MERVL-int flanked by two MT2_Mm LTRs) and solo LTRs (MT2_Mm not located in proximity offull-length MERVLs) in WT and Exosc3 cKO mESCs.

Figure 2.. Reduced Exosc3 levels in mESCs increase their developmental potential

(A) Volcano plot of differentially expressed genes and TEs between WT and Exosc3 cKO mESCs. 2CLC-specific genes and TEs, representing markers of 2CLCs, are highlighted in purple (defined from Eckersley-Maslin et al.). Other genes and TEs are marked in gray. (B) Gene set enrichment analysis (GSEA) of 2CLC genes/TEs in Exosc3 cKO mESCs. Genes/TEs are ranked according to the differential expression statistic (DESeq2 Wald test; lower rank, higher expression in cKO; higher rank, lower expression in cKO). Color bar displays differential expression statistic values. (NES, normalized enrichment score). (C) Uniform manifold approximation and projection (UMAP) (integrated data) of mESCs (n = 13,349 cells). Points (cells) are colored by annotated cell state. (D) Frequency of cell states in WT and Exosc3 cKO mESCs per biological replicate (*p < 0.05, **p < 0.005, Student’s t test; n = 3 per condition). (E) UMAP (integrated data) of mESCs. Points are colored by pseudotime (Slingshot trajectory analysis) and overlaid with RNA velocity vectors. (F) Distribution of cells along pseudotime from WT and Exosc3 cKO mESCs, grouped by biological replicate (differences in distributions were evaluated per replicate pair by Kolmogorov-Smirnov test; representative p value indicated from replicate 1). (G) Heatmap displaying H3K27ac pileup at enhancers in mESCs. (H) TF motifs significantly enriched in the mESC-specific and exosome-sensitive enhancer clusters from (from Figure 2G). Heatmap displays the −log₁₀(adjusted p value) of each target motif vs. background. (I) Aggregate signal at differential HiChIP loops (2D histogram). (J) Log₂ enrichments (observed/expected) for the overlap of HiChIP loop anchors and enhancer clusters or 2CLC gene TSSs (***p < 0.001). (K) Correlation of differential gene and TE expression in Exosc3 cKO and siExosc3 mESCs. (L) Schematic of experimental design for embryology experiments. (M) Bar plot displaying the number of dissected embryos containing embryonic only or embryonic and extraembryonic contribution from control siRNA and Exosc3 KD mESCs, respectively. (N) Embryo injected with Exosc3-siRNAs-transfected cells showing contribution to embryonic (epiblast) and extraembryonic (VYS, visceral yolk sac) compartments. White-dotted line outlines embryonic tissues at E6.5 not covered by extraembryonic membranes. (O) Embryos injected with Exosc3-siRNAs-transfected cells showing contribution to embryonic (NE, neuroectoderm) and extraembryonic (VYS, visceral yolk sac) compartments. Gray-dotted line outlines the silhouette of embryonic tissues at E8.0 covered by extraembryonic membranes. Bottom panels are at higher magnifications.

Figure 3.. Loss of Exosc3, Zcchc8, Rbm7, and Ints11 causes mESCs to acquire a 2CLC gene signature

(A) Analysis workflow for the integrated analysis of bulk RNA-seq and TT-seq data. (B) Density plots displaying the distribution of log₂-fold changes from a differential expression analysis of TT-seq labeled vs. unlabeled RNA in WT and Exosc3 cKO mESCs. Positive-fold changes indicate higher expression in labeled RNA. (C) Density plots displaying the distribution of log₂-fold changes of genes from a differential expression analysis of TT-seq labeled vs. unlabeled RNA in WT mESCs, grouped by the number of exons. Positive-fold changes indicate higher expression in labeled RNA. (D–F) GSEA of 2CLC genes/TEs in differential expression signatures from Zcchc8 KO, Rbm7 KD, and Ints11 KD mESCs.

Figure 4.. Gene expression asymmetry upon loss of Exosc3, Zcchc8, Rbm7, and Ints11

(A and B) Density plots displaying the differential expression statistic (Wald test) from a differential expression analysis of Exosc3 cKO, Zcchc8 KO, Rbm7 KD, and Ints11 KD mESCs, grouped by gene biotype. Protein-coding genes are further grouped by the number of exons. Fill colors represent median differential expression statistic values per group. (C) Scatter plot displaying correlation of differential expression statistic values between WT vs. Exosc3 cKO (x axis) and Zcchc8, Rbm7, and Ints11 KD gene expression signatures (y axis) in mESCs. 2CLC genes/TEs are displayed in purple, and other genes/TEs are displayed in gray. (D) GSEA of 2CLC genes/TEs in differential expression signatures from Exosc3 rapid degradation (FKBP) in mESCs at 6 and 12 h. (E) Box plots displaying the differential expression statistic from a differential expression analysis comparing 6 and 12 h Exosc3-depleted mESCs for selected gene biotypes. p-values were calculated using an unpaired, two-sided Wilcoxon rank-sum test and adjusted using the Benjamini-Hochberg method (*p < 0.05, ****p < 0.0001).

Figure 5.. The nuclear RNA exosome controls premature termination

(A) Pausing ratio of RNAPII ChIP-seq in WT and Exosc3 cKO mESCs. (B) As above, but with metabolically labeled RNA. (C and D) Metagene plots of RNAPII ChIP-seq with 8WG16 (top) and CTD4H8 (bottom) antibodies in WT and Exosc3 cKO mESCs, grouped by exon count. (E) Boxplots displaying normalized levels of labeled RNA in gene bodies in WT and Exosc3 cKO mESCs, grouped by exon count. p-values were calculated using an unpaired, two-sided Wilcoxon rank-sum test and adjusted using the Benjamini-Hochberg method (*p < 0.05, ****p < 0.0001).

(F) Metagene plot displaying the coverage Zcchc8 eCLIP peaks relative to the transcription start site (TSS) and transcription end site (TES). (G) Metagene plot displaying the frequency of unique pA+ and pA+,− termination sites in WT mESCs relative to gene TSS and TESs. (H) GSEA of genes with premature termination events in WT mESCs, calculated in differential expression signatures derived from Exosc3 cKO mESCs. (I) As above, but signatures derived from Exosc3-depleted mESCs by FKBP at 6 (center) and 12 h (right). (J) Bar plot displaying normalized enrichment scores (NESs) and p values (adjusted using the Benjamini-Hochberg method) from GSEA of genes with premature termination events in WT mESCs across differential expression signatures. (K) GSEA of 2C-specific genes and blastocyst genes with and without premature termination events in a differential expression signature from WT vs. Exosc3 cKO mESCs.