The corepressor NCOR1 and OCT4 facilitate early reprogramming by suppressing fibroblast gene expression

Data availability

New RNA-sequencing data generated for this article can be accessed via NCBI GEO database with the accession number GSE139376. High-content screening data can be downloaded as supplemental data from our previous study (Peñalosa-Ruiz et al., 2019). Other RNA-sequencing data generated in our lab previously, include the reprogramming siRNA non-targeting controls and time-course RNA-sequencing in reprogramming (Peñalosa-Ruiz et al., 2019). These can be accessed under a superseries with accession number GSE118680.

MEF-to-iPSC reprogramming

Primary mouse embryonic fibroblasts (MEFs) were derived from mouse wild type strain C57/BL6 (Erasmus MC iPSC facility). Passage 0–1 MEFs were seeded at a density of 10,000 cells per cm² in MEF medium consisting of 15% FBS, 100 nM β-mercaptoethanol and 1% non-essential aminoacids (Thermo Scientific, Waltham, MA, USA). Next day, MEFs were transduced with the doxycycline-inducible mouse tet-STEMCCA (Sommer et al., 2009) and rtTA lentiviruses (Addgene # 20342) at an MOI of 1. One day after transduction, cells were either transfected with siRNAs or induced for reprogramming with reprogramming medium. This consisted of DMEM-high glucose (Thermo Scientific, Waltham, MA, USA) supplemented with 10% stem cell grade FBS (Hyclone, Logan, UT, USA), 200 nM β-mercaptoethanol (Sigma, Kawasaki, Kanagawa, Japan), 1% sodium pyruvate (Thermo Scientific, Waltham, MA, USA), 2 μg·mL⁻¹ doxycycline, 3 μM Chiron (GSK3-inhibitor), 0.25 μM Alk5i (TGF-ß inhibitor), 50 μg·mL⁻¹ ascorbic acid and 1 × 10³ U·mL⁻¹ LIF. Next day after adding reprogramming medium was considered reprogramming Day 1.

siRNA transfections

Transfections were performed as described previously (Peñalosa-Ruiz et al., 2019). Briefly, before adding the cell suspension containing OSKM-transduced MEFs, multi wells were prepared with transfection mix. This consisted of a 40 nM pool of 3 different siRNAs per target mRNA, diluted in Optimem (Thermo Scientific, Waltham, MA, USA) together with RNAiMAX lipofectamine, according to the manufacturers instructions. After this incubation, cell suspensions were added to each well. A total of 24 h after transfections, reprogramming started by adding medium with doxycycline.

High-content screening

High-content screening was performed as described before (Peñalosa-Ruiz et al., 2019). Briefly, transfected cells were cultured in Cell Carrier 96-well black plates (Perkin Elmer, Waltham, MA, USA), after 6 days, samples were fixed with 4% PFA and stained for CDH1 (Cell Signaling 14472) and SALL4 (29112; Abcam, Cambridge, UK) with DAPI counterstain. Plates were imaged in an Opera High-Content-Screening System (Perkin Elmer, Waltham, MA, USA). Colony segmentation was done in multiple Z-planes using SALL4 staining. Columbus Software (Perkin Elmer, Waltham, MA, USA) was used to extract all features in an automated fashion. Z-score normalisation per plate was applied. The normalised data for high-content features were used for knockdown-to-knockdown Pearson-correlation analysis.

RNA isolation and RT-qPCR

RNA was isolated with RNA Micro-prep kit (Zymo Research, Irvine, CA, USA). Integrity was verified with Bioanalyzer and concentrations were determined with Nanodrop. For reverse transcription we used SuperScript III Kit (Thermo Scientific, Waltham, MA, USA) starting with 120–180 ng total RNA. For the RT-qPCR reaction 1–2 ng of cDNA were used in 20 μL reaction mix with Sybr Green mix ready to use (iQ-SYBR-Green Supermix; Biorad, Hercules, CA, USA). Relative gene expression was calculated with the ΔΔCt method using GAPDH as reference.

RNA-sequencing library preparation and data analysis

Kapa-RNA HyperPrep kit with RiboErase (Roche; Kapa Biosystems, Wilmington, MA, USA) was used for ribosomal RNA depletion and library preparation, with 200 ng total RNA. Library amplification was performed for 10 cycles, according to manufacturer’s instructions. Correct size distribution of 300 bp was verified with Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Quality controls were performed for selected transcripts in each sample before and after library preparation by qPCR.

NextSeq500 Illumina platform was used for paired-end library sequencing, with 43 bp read length. STAR v.2.5.b (Dobin et al., 2013) was used to align reads to mouse genome assembly mm10. Data was normalised to log2-cpm with R package edgeR v.3.20.9 (Robinson, McCarthy & Smyth, 2010).

For differential gene expression analysis we used R Package DESeq2 v.1.18.1 (Love, Huber & Anders, 2014). Gene lists were filtered with a cut-off of p-adjusted value < 0.05.

Gene Ontology classification with differential gene lists was performed with web-based DAVID (Huang, Sherman & Lempicki, 2009), considering only categories with gene counts > 10 and p-adjusted value < 0.05.

ChIP-seq and RNA-seq data integration

For ChIP-seq data integration, FASTQ files of ChIP-seq data were downloaded from NCBI GEO (GSE70736, GSE90893, GSE90895) and mapped to the mm10 genome assembly using BWA allowing one mismatch per read. NCOR1/NCOR2 peaks (Zhuang et al., 2018) were intersected with siNcor1 day 6-upregulated genes (linux join, using gene symbols) to obtain genomic locations of interest. Heatmaps of ChIP-seq data were generated with fluff (Georgiou & Van Heeringen, 2016) with hierarchical clustering and RPKM normalisation. RNA-seq data (Chronis et al., 2017) was intersected with differentially expressed genes in our study (linux join, using gene symbols). Corresponding RNA-seq heatmaps were generated using seaborn clustermap with row clustering and row z-score normalisation. Venn diagrams were obtained using venn3 from the matplotlib_venn library. Hypergeometric probabilities were calculated using scipy.stats.

High-content screening reveals early reprogramming disruption upon Ncor1 depletion

In a previous study, we performed a high-content siRNA screen, targeting 300 chromatin-associated factors during early reprogramming (Fig. 1A) (Peñalosa-Ruiz et al., 2019). We hypothesised that the onset of reprogramming would be associated with colony-level phenotypic changes and that these changes are affected upon gene perturbation. Moreover, using computational analysis of high-content imaging data, perturbed genes can be grouped according to their phenotypic similarities in a multidimensional space, rather than focusing on a single phenotype (Boutros, Brás & Huber, 2006). We used early pluripotency markers CDH1 and SALL4 to detect colonies at day 6 of reprogramming. Additionally, multiple colony features were extracted, including pluripotency marker expression levels, morphological traits such as colony roundness, area and symmetry, and many other colony texture and morphological features (Peñalosa-Ruiz et al., 2019). We selected 30 hits that were subjected to a secondary RNA-seq-based screening (Peñalosa-Ruiz et al., 2019) (overview in Fig. 1A). We have previously identified functional interactions from knockdowns showing similar phenotypes (Mulder et al., 2012; Peñalosa-Ruiz et al., 2019; Tanis et al., 2018). Using this rationale, we compared knockdown-to-knockdown correlations based on high-content phenotypes and also based on transcriptomes (Fig. 1B). We filtered for pairwise correlations showing highest scores in both the high-content and the RNA-seq screen (both R > 0.4). We visualised the resulting potential interactions between genes (nodes) as networks, considering both high-content imaging correlations (edge width) and transcriptome-based correlations (edge colour) of all siRNA knockdowns. We previously uncovered a functional interaction between BRCA1, BARD1 and WDR5 (Peñalosa-Ruiz et al., 2019). The strongest correlation however is between NCOR1-OCT4, raising the possibility of a functional interaction between these factors (Fig. 1B).

We verified that, based on high-content features, siNcor1 positively correlated with siOct4, but anti-correlated with the negative (nt, non-targeting) control (Figs. 1C and 1D). Colony number analysis derived from the screen showed a significant reduction of colony formation upon Ncor1 depletion (Fig. 1E; Fig. S1). High-content images from Ncor1 and Oct4 knockdowns showed very low CDH1 and SALL4 intensities and compromised colony formation, when compared to the nt control. We validated this finding using an independent In-Cell Western assay and confirmed that siNcor1 and siOct4 show similarly impaired reprogramming phenotypes (Fig. 1F; Fig. S1).

From these experiments we concluded that siRNA targeting of Ncor1 results in an early reprogramming phenotype that is similar to that caused by Oct4 knockdown.

Ncor1 depletion affects transcriptional regulation of fibroblast identity

To gain insight in the molecular events that explain early reprogramming disruption by Ncor1 depletion, we performed gene expression analyses. We first asked whether Ncor1 mRNA expression followed a dynamic pattern during early reprogramming. For this, Ncor1 mRNA expression was quantified by RT-qPCR at different reprogramming timepoints, from MEFs to reprogramming day 6. We found that Ncor1 mRNA expression was relatively stable throughout reprogramming, with slight variations (Fig. 2A). Then, we sought to understand the role of NCOR1 in transcriptional dynamics during early phases of iPSC reprogramming. For this purpose, we performed RNA-sequencing at day 3 and day 6 after expression of OSKM factors in MEFs transfected with control or Ncor1 targeting siRNAs, respectively. Differential gene expression analysis revealed 405 deregulated genes at day 6 (nt control vs siNcor1, adjusted p-value < 0.05, Figs. 2B and 2C). The effect size was relatively small for most transcripts, suggesting that the strong defects in reprogramming observed with siNcor1, are brought about by moderate changes in expression of hundreds of transcripts (mean absolute log₂-fold change of 0.43). The majority of the genes (94%) were up-regulated, as expected for a transcriptional corepressor, although the deregulation of some of these genes may be an indirect consequence of the reduction of NCOR1 rather than a reduced NCOR1-mediated repression on that particular gene.

To gain insight into which processes were affected in Ncor1-depleted cells, we performed Gene Ontology classification (GO) for upregulated genes because the downregulated were very few (Fig. 2D). This analysis showed that upregulated genes were associated with cell–cell adhesion, poly-A-RNA-binding, exosome and cytoskeleton (Fig. 2D). Using our reprogramming RNA-seq time course dataset (Peñalosa-Ruiz et al., 2019), we observed that genes associated with these terms were indeed downregulated in the course of early reprogramming (Fig. 2E). To note, the cell adhesion terms featured laminins, collagens and other genes related to basement membrane and extracellular matrix. These data suggest that reduction of Ncor1 prevents downregulation of genes related to fibroblast identity during MEF to iPSC reprogramming.

NCOR1 and OCT4 co-regulate two distinct MEF-expressed groups of genes

To gain insight into the functional relationship between OCT4 and NCOR1, we examined RNA-seq profiles at reprogramming day 6 for both knockdowns. To verify Ncor1 and Oct4 mRNA-targeting by the corresponding siRNAs, we plotted the RNA-seq normalised values for each of the knockdowns, compared to the controls at day 3 (Fig. S2). Oct4 knockdown did not show a high efficiency of knockdown at day 3 in these samples, however from time course experiments we have observed that Oct4 mRNA is targeted and efficiently silenced at day 2, but rapidly recovers within a couple of days (Fig. S2). This is due to its high overexpression driven by the tet-OSKM cassette (Sommer et al., 2009). Even with such transient early knockdown we observe a strong phenotype (Fig. 1) and many deregulated genes (Fig. 2A), highlighting the important contribution of OCT4 in OSKM reprogramming.

When we compared genes deregulated in each knockdown we found that 43% of the genes deregulated in siNcor1 cells were also deregulated in siOct4 cells (3.4 fold enrichment compared to random overlap, hypergeometric p-value 8.9 × 10⁻⁵⁴; Fig. 3A). This corresponds to only ~5% of genes differentially expressed upon Oct4 knockdown. After hierarchical clustering of genes deregulated by both siNcor1 and siOct4 based on RNA-seq, we identified three clusters (cluster 1, 2A, 2B; Fig. 3B). Cluster 1 genes were relatively unaffected at day 3, but at day 6 they were downregulated for both knockdowns. For cluster 2A we found that genes were also unaffected at day 3, but showed a similar upregulated pattern in both Oct4 and Ncor1 siRNAs at day 6 (Fig. 3B). Cluster 2B contained genes with high expression in MEFs, most of which failed to be downregulated at day 3, whereas by day 6 they predominantly showed opposite expression changes in the knockdowns (siNcor1 vs siOct4) (Fig. 3B). Our RNA-seq time course dataset (Peñalosa-Ruiz et al., 2019) showed that genes in cluster 1 mostly remain the same in time, whereas clusters 2A and 2B tend to go down in time, but with different dynamics (Fig. 3C).

To shed some light on the difference between clusters 2A and 2B, we performed GO classification for genes of each cluster (Fig. 3D). NCOR1-OCT4 co-regulated genes (cluster 2A) showed an association with RNA metabolism and cell adhesion (Fig. 3D). The term “poly(A)-RNA binding” contained several components for pre-mRNA splicing. Genes increased by siNcor1 but not siOct4 (cluster 2B), were associated with functions in extracellular matrix organisation and other cell adhesion features. We wondered whether the “cell adhesion” terms were related to the mesenchymal-to-epithelial transition (MET) observed in early reprogramming (Li et al., 2010; Samavarchi-Tehrani et al., 2010). We found that most mesenchymal genes exhibit a cluster 2B-type pattern, with increased expression at day 3 in both knockdowns and opposite effects of the two knockdowns at day 6. Epithelial genes on the other hand show more variable patterns (Fig. 3E). In siOct4 cells, the downregulation of mesenchymal gene expression is delayed. By day 6, it not only catches up, this decrease in expression is more severe in siOct4 cells relatively to control cells, suggesting that OCT4 moderates this decrease during normal reprogramming. In siNcor1 cells on the other hand, mesenchymal gene expression is maintained at higher levels at both time points, although some epithelial genes are also expressed at higher levels (Fig. 3E).

These results identify two distinct groups of somatic genes, the expression of which is normally reduced during reprogramming. Ncor1 depletion prevents downregulation of both groups, whereas Oct4 depletion prevents the downregulation of one group and causes a stronger decrease in expression in the other group of genes.

NCOR1 and OCT4 target promoters associated with cellular and metabolic processes in MEFs and early reprogramming cells

To gain insight into the functional relationship between NCOR1 and OCT4, we analysed published ChIP-seq and RNA-seq data (Chronis et al., 2017; Zhuang et al., 2018). First, we asked whether the up-regulated genes that we identified in siNcor1 and siOct4 cells, were either MEF genes or genes dynamically changing during reprogramming. Consistent with our earlier analysis (Fig. 3C), these independent data confirmed that upregulated genes in siNcor1 cells were mostly MEF-expressed and very early reprogramming (48 h) genes, but were mostly inactive in ES cells (Fig. 4A). In the heatmap several genes related to fibroblast identity are labelled (lamins, collagens, fibronectin, signalling genes; Fig. 4A). In contrast, genes upregulated in siOct4 cells showed a variety of expression patterns during reprogramming (Fig. 4B).

Then, we asked which fraction of differentially expressed genes in Ncor1- or Oct4-knockdowns were bound by NCOR1 protein. We overlapped NCOR1 ChIP binding sites in reprogramming cells with siNcor1 and siOct4 differential genes (Figs. 4C and 4D) and calculated the fold enrichment. We found that NCOR1 binds to a significant fraction of deregulated genes in Ncor1 and Oct4 knockdowns (Figs. 4C and 4D; Table 1). It should be noted that the number of genes identified may be affected by the use of siRNA knockdown rather than genetic knockouts and the different time-points and different reprogramming conditions.

To gain insight on the function of transcripts affected by Ncor1-depletion in reprogramming, we exclusively selected NCOR1-bound genes that also showed increased expression after Ncor1 siRNA knock-down (see intersection of 121 genes in Fig. 4C). We used publicly available ChIP-seq datasets to ask whether OCT4 and c-MYC binding was dynamic at these genes during reprogramming (Figs. 4E–4G; Fig. S3). Quantifying the ChIP-seq data (reads per million) showed that c-MYC and OCT4 signal at these 121 genes is higher at 48 h after the onset of reprogramming than it is in mouse embryonic fibroblasts (MEFs) and Embryonic Stem Cells (ESC) (Fig. 4E, MEF data for OCT4 not available (N/A)). This indicates that during early stages of reprogramming OCT4 and c-MYC show increased binding to NCOR1 targets affected by Ncor1 knockdown. Interestingly, this binding activity is not sustained through the pluripotent state (ESC), hinting at a transient role for NCOR1 in reprogramming.

To investigate this process further, we analysed several chromatin features for the genes in question. In general, these genes contain the active gene marks H3K4me3 (Fig. S3A), indicating some transcriptional activity. Notably, we observed a gradual decrease in signal of the H3K27ac active modification when comparing MEFs, 48 h of OSKM-induced reprogramming and ESCs (Fig. 4G). Moreover, this decrease in H3K27ac levels was concurrent with reduced occupancy of the P300 histone acetyltransferase on these NCOR1 regulated genes (Fig. 4F middle panel). Conversely, binding of histone deacetylase HDAC1 is highest at 48 h after OSKM induction, compared to either MEFs or ESCs (Fig. 4G). These data support the notion that NCOR1 regulated genes are transcriptionally downregulated during the early phases of reprogramming. Additionally, the underlying mechanism may include NCOR1-mediated recruitment of HDAC proteins, including HDAC1. Possibly, reprogramming factors OCT4 and c-MYC could be assisting NCOR1 binding to such genomic targets. However, more experimental evidence is needed to test this.

To explore the function of NCOR1 targets also bound by OCT4 and c-MYC, we performed a Gene Ontology classification of genes analysed in Fig. 4D (Fig. S3). This analysis confirmed once more the terms we observed before (Fig. 2D; Fig. 3D), related to cell-cell adhesion, extracellular exosome, extracellular matrix, cytoskeleton, focal adhesion and poly(A) RNA binding, among others (Fig. S3). As before, the poly(A) RNA binding term included several genes related to translation and pre-mRNA splicing.