We used a machine learning approach, gapped-kmer SVM (gkm-SVM) (Ghandi et al., 2014a; Beer et al., 2020; Lee et al., 2015) to identify TFBS driving variations in chromatin accessibility during cell state transitions. We chose gkm-SVM for its combination of predictive accuracy and interpretability (Ghandi et al., 2014a; Beer et al., 2020; Lee et al., 2015; Shigaki et al., 2019; Beer, 2017; Yan et al., 2021; Kreimer et al., 2017; Gate et al., 2018; Mo et al., 2016; Jain et al., 2024), and because of its advantages relative to deep neural networks (DNNs) when training on smaller numbers of differentially active peaks during a transition between related cell types (∼2000 peaks). However, we have shown that gkm-SVM and DNNs detect similar features when DNNs produce a robust model (Shigaki, 2024). We trained a gkm-SVM model on the top 10,000 distal ATAC-seq peaks against random negative genomic sequences, following our standard pipeline (Beer et al., 2020), for ESC and DE (after 2 days of differentiation, DE-D2). Each gkm-SVM model summarizes the TFBS required to classify the ATAC-seq peaks in a vector of gapped kmer weights which can predict the accessibility of the sequences. These weight vectors have a tail of large positive weight kmers which map to activating TFs expressed or active in the cell type, and a much weaker tail of negative weight or repressive TFs, as shown in Figure 2A for ESC and DE. In Figure 2A, we also show the average of all weight vectors trained from every ENCODE DNase-seq dataset (1270 cell/tissues), and again find an imbalance between activating and repressive TFBS. Within the repressive tail of negative weight kmers learned in the full ENCODE set of tissues, we do learn kmers mapping to known repressive binding sites, most commonly SNAI1/2, TWIST, and ZEB1/2 (5′-CACCT-3′ motif) which are repressors known to play a key role in EMT and cancer (Craene and Berx, 2013; Vandewalle et al., 2008); GLI3, GLIS2, and ZBTB7A (5′-GACCCC-3′ motif) known to have a repressive role in Hedgehog signaling (Lex et al., 2022; 2020) and fetal hemoglobin switching (Liu et al., 2021); and nuclear hormone receptor monomer and dimer elements (5′-AGGTCA-3′) bound by large families of co-repressors (Perissi et al., 2010). We also learn non-ZEB-like E-box motifs (5′-CAGNTG-3′) which are usually activating, but repressive in heart and fibroblasts, where they could be bound by E-box binding helix-loop-helix TFs (e.g., E47, E2A, also known as TCF3) and targeted by repressive ID factors (Perk et al., 2005). Although we detect several well documented repressive TFBS, in all of our ENCODE models the negative weight tail of repressive TFBS is only slightly enriched above the expected null normal distribution, while the positive weight tail is dramatically enriched in all cell/tissues and contains hundreds of activating cell-specific TFs (5-10 per cell/tissue). In addition, we find the repressive TFBS in broader sets of cell types and tissues than the activators, so while they are clearly playing an important repressive role, there may not be cell-specific repressors unique to each cell type. For the specific case of the ESC to DE cell-state transition, we trained models on ESC and DE ATAC-seq peaks, and extracted motifs from each weight distribution (Methods). We find motifs for the expected regulators OCT4-SOX2-NANOG in ESC and GATA in DE (Figure 2B), but we found no or very weak negative weight motifs. Consistent with the results from our machine learning models, a systematic study of TF binding in HepG2 cells found that 34 out of 35 TFs with strongly reproducible effects on expression were activating, while only REST was repressive (fraction significant >0.8, |median effect estimate|>0.5) (Moyers et al., 2023). The miRNA repression analyzed in this paper is an alternative repressive mechanism.

We next trained a model on the most differentially active peaks in ESC vs. DE (n = 3000 in each set). Now OCT4-SOX2-NANOG is learned as a positive signal, and GATA, EOMES, SMAD, and MIXL1 (Crx homeodomain motif) are learned as negative signals. This training design (ESC peaks vs. DE peaks) more sensitively detects the motifs required to explain all the accessibility changes induced during the cell-state transition from ESC to DE. Again, the known activating TFBS signals OCT4-SOX2-NANOG and GATA, EOMES, SMAD, and MIXL1 are learned (the same motifs learned when trained against negative genomic regions are learned as contributing positively to that cell type’s differentially active peaks). If repressive TFs were contributing broadly to changes in ESC and DE chromatin accessibility, we would expect to see them. In some cases, we do learn an EBOX motif (5′-CANNTG-3′, which as discussed above can be repressive in some cell types) as contributing to ESC specific peaks, but this is most likely bound by MYC, which is strongly expressed in ESC but is downregulated in DE. This analysis does not preclude the possibility that DE TFs are activating a repressive TF that is only binding at a few core ESC TF enhancers (those flanking POU5F1 (OCT4), SOX2, and NANOG), but if this were a dominant mechanism, we would also expect to see the repressive TF having an effect at other peaks which have weak binding sites for the TF, and thus should appear in the weight vector. In conclusion, we see little evidence for a strong repressive TF expressed in DE which shuts down the ESC core TF enhancers, and little evidence for a strong repressive TF expressed in ESC which keeps the DE TFs repressed. The PRC2 complex may play this repressive role by being recruited to specific promoters, but how the complex achieves cell-specificity is unclear (Youmans et al., 2021).

The gene regulatory network model (Luo et al., 2023) we developed to understand the dynamic effect of CRISPRi enhancer perturbation is a two-state transcriptional model with a large enhancer activated transcription rate, e 1 , and a low basal transcription rate, e 0 (see Figure 3A for GRN model variables). When we learn the core TFs in any cell type with gkm-SVM, we also find binding sites for multiple core TFs in each core TF enhancer. For simplicity, we model this cooperative autoregulatory activity by assuming that all core TFs are activated together, and we call this activation ψ 1 t for cell type 1 and ψ 2 t for cell type 2. ψ t describes how the concentration of core TFs for each cell type vary as a function of time (Methods). Each cell type’s enhancers are active proportional to c ψ 3 b + c ψ 3 and inactive proportional to b b + c ψ 3 . We will assume that miRNAs are similarly transcriptionally activated by enhancers, with the same rates, and repress the other cell type’s core TF activity with an additional rate d . Relaxing the assumption that core TFs and miRNAs are transcribed at the same rates e 1 and e 0 does not qualitatively change our findings. The translational repression rate, d , includes miRNA processing and nuclear transport. Together this yields two coupled nonlinear differential equations for ψ 1 and ψ 2 as shown in Figure 3A. We have also added a stimulus, δ t , which acts at a subset of enhancers to induce the transition from cell type 1 to 2 by activating enhancers of cell type 2 through a mechanism independent of ψ 2 . Further details of the model are provided in Methods.

In the limit d = 0 , these two equations are completely uncoupled, and each has a stable OFF state at ψ = e 0 / r and a stable ON state at ψ = e 1 / r . When initial conditions are near the OFF state, a stimulus is required to transition into the ON state, but for d = 0 , and even for small d , it is possible to activate ψ 1 and ψ 2 at the same time. As discussed in the introduction, this likely causes inappropriate expression of many downstream genes and interferes with normal cell function and viability. For small d = 0.1 , with a stimulus δ 0 = 0.3 (see Methods) sufficient to induce the activation of cell type 2, ψ 2 becomes activated (Figure 3B, left), but cell type 1 remains active. In contrast, for larger miRNA repression of cell type 1’s TFs d = 1.0 , at the same stimulus ( δ 0 = 0.3 ), when ψ 2 becomes activated, the activation of miRNAs repressing ψ 1 TFs are sufficient to drive ψ 1 fully to the OFF state (Figure 3B, right).

Further understanding of this behavior is provided by phase plane analysis in Figure 3C. For d = 0.1 , there is a stable fixed point at ψ 1 , ψ 2 = 2.3,2.3) (Figure 3C, left), showing that a stable solution exists with both ESC and DE TFs activated. With sufficiently strong negative feedback, this stable fixed point disappears ( d > 0.4 ), and the only three stable fixed points are with ψ 1 in the ON state and ψ 2 in the OFF state ψ 1 , ψ 2 = 3,0), ψ 1 in the OFF state and ψ 2 in the ON state ψ 1 , ψ 2 = 0,3), and both ψ 1 and ψ 2 in the OFF state ψ 1 , ψ 2 = 0,0) as shown in Figure 3C (right), for d = 1 . This situation is more similar to what is seen experimentally in the ESC to DE transition, where activation of core DE TFs genes (GATA6, MIXL1, EOMES, SOX17) induces the inactivation of the core ESC TF genes (POU5F1, SOX2, and NANOG).

Next, we analyze miRNA expression data to detect miRNAs that are differentially expressed in ESC and DE and might be repressing TFs in the other cell type. PCA of miRNA expression across 177 ENCODE samples (Reese et al., 2023) shows some degree of clustering and separation between samples originating from the three germ layers (endoderm, mesoderm and ectoderm) (Figure 4A; Supplementary Figure S2, Methods). The three human embryonic stem cell samples (hESC) are close to each other in the PCA space, and have an average expression correlation of 0.95 compared to average expression correlation of 0.77 for all 177 pairwise comparisons.

We performed a differential expression analysis for miRNAs in DE compared to the ESC cell line and found 59 miRNAs overexpressed in DE state and 100 miRNAs with higher expression in ESC (Figures 4B, C). Some miRNAs upregulated in either ESC or DE are co-expressed across other tissues and cluster together (Figure 4B), providing evidence that miRNA expression differences are informative and vary across differentiation states and tissues.

ESC and DE ATAC-seq along with H3K27ac (Li et al., 2019; Luo et al., 2023) shows that peaks flanking differentially expressed miRNAs are consistently more active in the tissue where the flanking miRNAs are overexpressed (Figure 4D, Methods). Core DE and core ESC TFs bind preferentially to the most differentially accessible peaks between DE and ESC that are near a differentially expressed miRNA (Figure 4E). GATA6, EOMES and SOX17 have a significantly higher binding to ATAC peaks flanking miRNAs upregulated in DE (named “DE-MIR-flanking” peaks), whereas MYC, POU5F1 and NANOG have a higher ChIP-seq signal at differentially accessible peaks in ESC flanking ESC miRNAs (“ESC-MIR-flanking” peaks) (t-test p < 0.01, Figure 4E). In addition, GATA and POU5F1 DNA-binding motif sequences have a significantly higher enrichment in differentially accessible DE-MIR-flanking peaks and ESC-MIR-flanking peaks, respectively (t-test p < 0.05, Supplementary Figure S3, Methods). Taken together, there are distinct patterns of accessibility and TF binding at peaks near differentially expressed miRNAs in ESC and DE which suggest that they are functioning as enhancers driving the differential expression of these miRNAs in ESC and DE, consistent with our modeling in Figure 3.

Our dynamic network modeling (Figures 1, 3) highlights the possibility of miRNAs repressing core ESC and DE TFs during embryonic differentiation. We used TargetScan and mirDIP (McGeary et al., 2019; Tokar et al., 2018) to find the set of candidate miRNAs targeting 3′UTRs of core TFs (Figure 5A, Methods) and discovered a subset of miRNAs predicted to target MYC, SOX2, POU5F1, EOMES and GATA6 (Figure 5B). MIR92B is expressed at a higher level in ESC and is confidently predicted to target EOMES and GATA6 (Figure 5C). MIR145 is upregulated in DE and consistent with the target prediction, MIR145 was shown to repress pluripotency by targeting ESC TFs POU5F1, SOX2 and KLF4 (Xu et al., 2009). MIR892A is upregulated in DE while it has low expression in the ESC state and is predicted to target POU5F1. Another core ESC TF is MYC predicted to be targeted by MIRLET7D (Supplementary Figure S4). In addition, our analysis of differentially expressed miRNAs identified a potential combinatorial regulation of SOX2 by multiple DE miRNAs such as MIR145, MIR369 and MIRLET7F1 (Figure 5D). This potential combinatorial repression is consistent with SOX2 being the core ESC TF having the highest and the fastest downregulation following ESC-DE transition through DE day 3 (Supplementary Figure S1).

We identified differentially expressed miRNAs putatively targeting DE or ESC TFs and found differentially active enhancer peaks consistent with their cell-specific expression. MIR92B is predicted to target DE TFs and has a differentially active candidate enhancer (E1) flanking it (Figure 6A), where ESC TFs such as POU5F1 and NANOG bind. Conversely, DE TFs such as EOMES, GATA6 and SOX17 bind to the differentially accessible candidate enhancers (E2 and E3) flanking MIR145 and MIR892A (Figures 6B, C). These candidate enhancers are all within loop extrusion model predicted CTCF loops (Xi and Beer, 2021) enclosing the miRNA, providing more support that they are responsible for the miRNA differential expression. Together, these results suggest DE and ESC TFs, through activating cis-regulatory elements, upregulate miRNAs, which in turn repress core TFs for the other transition state.

We trained gkm-SVM on ATAC-seq data from ESC and DE (Li et al., 2019), after removing promoter peaks within 2 kb of an annotated TSS, non-cell specific peaks accessible in more than 30% of ENCODE samples, and CTCF binding sites. We removed promoters as they have similar TFBS and are usually active across all cell types (Oh and Beer, 2024). Gkm-SVM uses counts of gapped kmer features under the 300bp sequences centered on the ATAC-seq peaks, as gapped kmers are effective at describing regulatory sequence features (Ghandi et al., 2014b; Amanchy et al., 2011) and do not require previous knowledge of TFBS. Following the default training method in the gkmSVM-R package (Ghandi et al., 2016), the AUROC was evaluated on 5-fold cross validation (CV) and AUROC is the average of 5 CV test sets (0.902 for ESC, 0.949 for DE-D2, and 0.975 for DE-D2 vs. ESC). All cross-fold validation sets produce similar sets of features and have very similar weight vectors. We normalized weight vectors to have zero mean and unit standard deviation for comparison. We extracted TFBS from this weight vector using gkm-PWM (Shigaki, 2024), which produces a lasso weight reflecting the importance of each motif in explaining the full weight vector (W in Figure 2). Other motif metrics Z and I describe the Z-score for kmers mapping to the motif and the error induced by removing the motif from the list.

We modified the gene regulatory network model (Luo et al., 2023) to include miRNA repression. The model describes the rate of change of cell-specific TFs, where one cell state’s TF concentrations are described by ψ 1 t , and the other cell state’s TF concentrations are described by ψ 2 t . For simplicity, we assume that through autoregulation, each cell state’s TFs turn on or off together, and can thus be represented with one concentration (Luo et al., 2023). The differential equation describes the change in concentration of the TFs, which is driven by the balance between a transcription rate, T ψ , and degradation/export, − r ψ , see, for example, (Alon, 2019). For both the TF and miRNA transcription rate we choose the form T ψ = e 1 c ψ 3 + e 0 b b + c ψ 3 . This rate depends on the TF activity, ψ , through enhancer activation, and the parameters control the rate at which transcription transitions from low transcription to high transcription as a function of ψ . It is a simple Hill form (Alon, 2019), but can also be understood as a two-state model with a large enhancer activated transcription rate, e 1 , and a low basal transcription rate, e 0 . Each cell type’s enhancers are active proportional to c ψ 3 b + c ψ 3 and inactive proportional to b b + c ψ 3 . We will assume that miRNAs are similarly activated with enhancers, and repress the other cell type’s core TF activity with a rate d . This rate includes miRNA processing and nuclear transport. The core TFs are degraded and/or exported from the nucleus at a rate r . Together this yields two coupled nonlinear differential equations for ψ 1 and ψ 2 as shown in Figure 3A. Transitions from cell type 1 to cell type 2 are modeled (Luo et al., 2023) by adding a time-dependent stimulus, δ t , which acts at a subset of cell type 2 enhancers to induce activation of ψ 2 , by increasing the probability that ψ 2 ’s enhancers are active. The probability that ψ 2 is in the basal state is b , and the probability that ψ 2 ’s enhancers are active is c ψ 2 3 without stimulus and c ψ 2 3 + δ t with the stimulus, with δ t = δ 0 1 + erf t − t 0 2 2 / 2 . This stimulus turns on gradually at t 0 = 3 . The other parameters are b , c , e 0 , e 1 , r , d , δ 0 = 0.5 , 1 , 0.1 , 3 , 1 , 1 , 0 unless specified otherwise.

We downloaded 177 ENCODE miRNA-seq samples from the ENCODE portal (Reese et al., 2023) (Supplementary Table S1). Raw miRNA expression counts were mean-normalized across the 177 samples and log₂-transformed. 835 miRNAs of the annotated 1877 miRNAs had a log₂-normalized expression above 4 (∼16 normalized counts) in at least two samples. We used this subset of 835 miRNAs to plot a PCA and heatmap of miRNAs across different tissues. We generated a table assigning a likely germ layer of origin to miRNA-seq tissue and cell line experiments (Supplementary Table S1). PCA and heatmaps were generated in R using ‘prcomp’, ggplot2 and pheatmap. We used ENCSR588ZRK (H1 embryonic stem cells) labeled “ESC” and ENCSR820BMF (endodermal cell originated from H1) labeled “DE” for miRNA expression analysis.

We used a threshold of |log₂FC| > 2 to call differentially expressed miRNAs between ESC and DE. Ensembl bioMart V.110 and V.111 were used for miRNA gene name conversion and transcription start site annotation (Martin et al., 2023).

ATAC-seq, ChIP-seq (TFs and H3K27ac) fastq and bigwig files were obtained from GEO GSE213394 (Li et al., 2019). They were processed according to the pipeline described by Shigaki et al. (Shigaki et al., 2019). We compared ATAC-seq and ChIP-seq samples in ESC to DE-D2 (day 2), as some TF ChIP-seq data was only available in DE-D2 (and not in DE-D3). We used ATAC-seq and H3K27ac to identify differentially active candidate enhancers regulating differentially expressed miRNAs in DE and ESC. We extended the ATAC peak summit by ± 500bp or 150bp on each side (peak length = 1000bp for H3K27ac and 300bp for ATAC) and calculated the average H3K27ac signal across the extended region, as H3K27ac peaks typically flank the ATAC peak summit. Differentially active regions were called from the set of top 25,000 strongest ATAC peaks in DE and ESC which satisfy the following: i) ATAC |log₂FC| > 1; ii) H3K27ac |log₂FC| > 0.25; iii) fall within 100kbp of a differentially expressed miRNA. Peak color indicates whether they are close to a miRNA upregulated in DE (“DE-MIR-flanking” peak) or to a miRNA upregulated in ESC (“ESC-MIR-flanking” peak). Of the 101 DE-MIR-flanking and 183 ESC-MIR-flanking peaks, only seven peaks are shared among the two sets and this subset of peaks has a relatively similar ATAC signal in DE and ESC (not differentially active and not affecting the rest of the analysis).

GATA6, EOMES and SOX17 ChIP-seq samples in DE were obtained from GSE213394, while the bigwig files for MYC (ENCFF145JGY, Snyder Lab), POU5F1 (ENCFF106YHB, Myers Lab) and NANOG (ENCFF512EZC, Myers Lab) ChIP-seq in H1 (ESC) were downloaded from the ENCODE portal (Luo et al., 2020). UCSC bigWigAverageOverBed was used to calculate average ChIP-seq signal over DE-MIR-flanking and ESC-MIR-flanking peaks and the values were log₂-transformed and t-test was used.

For motif enrichment analysis, we used STORM (Schones et al., 2007) and the following PWMs: GATA (MEME “MOTIF GATA1-ext”), GSC (JASPAR “MA0648.1”), SOX (JASPAR “MA0143.4”), EOMES (MEME “MOTIF Eomes”), NANOG (HOMER “motif224”), POU5F1 (JASPAR “MA1115.1”) (Bailey et al., 2015; Fornes et al., 2020; Heinz et al., 2010). t-test was used. To remove CTCF binding sites we scored with ScanACE (Karnik and Beer, 2015) with a score threshold of 5.0.

We used mirDIP to find a set of candidate miRNAs targeting core DE/ESC TFs and selected pairs of miRNA and core TFs whose mirDIP score is greater than 0.20. TargetScan was used to validate mirDIP predictions and align 3′UTRs to miRNA seed regions.

We used WashU Epigenome Browser (Li et al., 2022) for Figure 6. For better visual clarity in the figure, we used RefGene annotations (genomic coordinates) for MIR145 and MIR92B, while we used GENCODE V39 for MIR892A coordinates.