Reconstructing epigenomic dynamics through a single-cell multi-epigenome data integration framework

Domain: neuroscience
Published: 2025-12-17
Editorial outcome: minor_revision

Abstract
Transcriptional regulation arises from the dynamic and combinatorial actions of multiple regulatory factors on genomic DNA. Although many epigenomic regulators have been identified, the precise order in which these factors accumulate at individual gene loci to activate transcription remains unclear. Here we show a single-cell data integration framework that infers the binding order of multiple chromatin factors at single-cell resolution. Central to this framework is sci-mtChIL-seq, a scalable single-cell method that simultaneously profiles genome-wide binding of RNA polymerase II (RNAPII) and diverse epigenomic regulators. By defining transcriptional states through RNAPII occupancy and integrating multiple sci-mtChIL-seq datasets, we systematically link the combinatorial patterns of transcription factor binding, histone modifications and chromatin remodeling. This framework reveals the temporal coordination among chromatin factors during transcriptional activation, providing a powerful approach to uncover context-dependent epigenomic dynamics and the principles of gene regulation in complex cellular systems.

Introduction
Changes of cell state induced by differentiation signals, ageing, or stress are fundamentally driven by changes in gene expression. Epigenomic states, including chromatin structure and higher-order genome organization, play a central role in regulating gene expression. Comprehensive profiling of genome-wide epigenomic landscapes is thus essential to elucidate the mechanisms of transcriptional regulation1,2. However, epigenomic states are determined by the coordinated actions of numerous factors, including histone modifications and DNA-binding proteins. Despite increasing knowledge of individual regulators, the cooperative dynamics by which these factors collectively control transcription remain poorly understood.

To uncover such regulatory mechanisms, it is necessary to adopt an integrative approach that captures transcriptional state changes while simultaneously assessing the patterns of accumulation of multiple epigenomic factors. Two primary strategies have emerged to address this challenge. The first involves technological advances enabling the simultaneous profiling of an increasing number of factors3–9. For instance, uCoTarget9 has demonstrated the concurrent analysis of up to five factors. However, as the number of targets increases, the read depth per target tends to decrease, and physical interference between co-localized factors poses a fundamental limitation, making detection increasingly difficult.

The second approach leverages computational integration to generate datasets that allow in silico evaluation of numerous factors simultaneously10–12. Paired-Tag10, for example, measures histone modifications and mRNA expression concurrently, enabling integration based on transcript levels to identify multimodal epigenomic states. While this method has been applied primarily to assess the dynamics of accumulation at candidate cis-regulatory elements (cCREs), it highlights the potential to evaluate complex interactions between transcriptional states and multiple regulatory factors.

In this study, we developed a data integration framework to infer chromatin structural dynamics based on RNA polymerase II (RNAPII) binding sites. Building upon our previously established mtChIL-seq method13–15, we introduced sci-ChIL-seq for single-cell epigenome analysis, and sci-mtChIL-seq for single-cell multi-epigenome profiling. These methods enabled simultaneous mapping of genome-wide RNAPII binding and the positioning of various epigenomic factors, thereby linking epigenomic states to transcriptional activity. By assuming that cells with similar transcriptional states exhibit comparable epigenomic landscapes16, we performed RNAPII-anchored data integration across multiple datasets. This allowed us to propose a framework that visualizes and infers the dynamic epigenomic regulation orchestrated by multiple factors during transcriptional activation.

Results
Development and validation of the sci-ChIL-seq method
To infer chromatin structural dynamics associated with transcriptional states, we developed an integrative framework that defines cell states based on RNAPII binding and enables the unified analysis of multiple epigenomic factors. A key requirement for this approach was the ability to simultaneously profile genome-wide binding of RNAPII along with other chromatin regulatory factors, such as histone modifications, at the single-cell level. Existing methods often lack the capacity to capture both transcriptional activity and diverse epigenomic landscapes concurrently in individual cells. Therefore, building on our previously established ChIL-seq platform13,14, we developed a single-cell multi-epigenome analysis method to overcome these limitations and enable comprehensive, simultaneous profiling of RNAPII and other key epigenomic determinants.

First, we combined the split-and-pool approach17 with ChIL-seq to develop a high-throughput sci-ChIL-seq for single-cell-level analysis of the binding state of a single factor (Fig. 1a). ChIL-seq uses ChIL probes (secondary antibodies conjugated with oligonucleotides) that contain Tn5 transposon target sequences and T7 promoter elements, allowing localized genomic regions near primary antibody binding sites to be tagged, transcribed, and amplified for sequencing.

In sci-ChIL-seq, cell barcoding is achieved via two rounds of indexing separated by a split-and-pool step17,18. sci-ChIL probes with unique index sequences are introduced and incorporated near target sites via Tn5-mediated transposition, assigning the first index. Prior to this, cells are fixed, permeabilized, and blocked following our previously established ChIL-seq protocols13,14. Cells are then detached from the plate, pooled, and redistributed into a second 96-well plate, where indexed genomic regions are transcribed and reverse-transcribed to cDNA, followed by a second indexing via PCR. We distributed 25 cells per well during split-pool barcoding to reduce the probability that two cells with the same T5 index are sorted into the same well. This approach minimizes index collisions and doublet events, and is a standard practice in sci-based methods to preserve single-cell resolution.

Using 96 unique sci-ChIL probes and redistributing 25 cells per well, up to 2400 single-cell datasets can be obtained per plate. The sci-ChIL probe sequence includes a T7 promoter, a partial universal adapter19,20, a 6 bp T5 index, extended mosaic ends, and Tn5 binding sequences, totaling 80 bp (Supplementary Fig. 1a).

By contrast, the sci-ChIL probe reaction must be performed with the cells fixed on the plate. To ensure the success of on-plate reactions, cells were pre-labeled with biotin-SS-sulfo-OSu, enabling firm attachment to NeutrAvidin-coated plates during the first indexing step. Cells were then efficiently detached using reducing agents and protease Accumax (Supplementary Fig. 1b–d), as validated with NIH3T3 fibroblasts. Biotinylated cells adhered stably without aggregation, whereas non-biotinylated cells were lost during washing. Detachment was efficiently achieved through chemical treatment.

We next applied sci-ChIL-seq to measure RNAPII binding and evaluate transcriptional states. Using NIH3T3 cells, we assessed data quality by determining the number of unique fragments per cell for RNAPII and H3K4me3 targets (Fig. 1b). Median values were 5163 and 5877 for RNAPII, and 7534 and 9039 for H3K4me3. These median values were comparable or superior to those from other methods, such as NTT-seq (382–2854 per target), MulTI-Tag (193–3568 per target), and uCoTarget (1390–11,437 per target).

We further tested whether RNAPII signals obtained from sci-ChIL-seq could serve as proxies for transcriptional states similar to RNA-seq. Bulk RNA-seq data from NIH3T3 cells were used to group genes into quintiles by expression level. The pseudo-bulk of the sci-ChIL-seq data showed that RNAPII and H3K4me3 accumulation at gene bodies or transcription start site (TSS) regions, respectively, increased with higher gene expression (Fig. 1c, Supplementary Fig. 1e)21, with RNAPII exhibiting a near-linear relationship. At the single-cell level, consistent binding patterns were observed across cells for both RNAPII and H3K4me3 within a defined genomic region (chr1:32.1–36.1 Mbp) (Fig. 1d). These results demonstrate that sci-ChIL-seq enables reproducible detection of factor-specific chromatin landscapes at single-cell resolution, and that RNAPII binding levels can be used to assess transcriptional activity in individual cells.

Simultaneous profiling by sci-mtChIL-seq
We next expanded sci-ChIL-seq to a multi-target format, developing sci-mtChIL-seq (Fig. 1e, Supplementary Fig. 1f), and implemented this using RNAPII and H3K4me3 in NIH3T3 cells. To confirm the fidelity of sci-mtChIL-seq, we compared the signals obtained with those from single-target sci-ChIL-seq. Fragment counts in 5 kb bins were correlated between pseudo-bulk datasets (Fig. 1f), showing high correlation coefficients (>0.8) between replicates and across both single- and multi-target experiments. Single-cell reproducibility was further assessed, with RNAPII and H3K4me3 signals showing factor-specific patterns similar to those observed in Fig. 1d, confirming that multi-target data matched the single-target results (Fig. 1g). Signal specificity was validated by comparing sci-mtChIL-seq datasets obtained with and without primary antibodies. Unique fragment counts per cell were higher in antibody-treated samples, indicating antibody-dependent signal enrichment (Fig. 1h). Although fragment counts were reduced compared with those in single-target sci-ChIL-seq, they remained sufficient relative to other methods. To visualize signal specificity at the single-cell level, unique fragments per 5 kb bin were plotted using UMAP (Fig. 1i). Cells formed distinct clusters based on the presence or absence of primary antibodies, with matched data from the same cell linked and color-coded. This demonstrated a clear separation of specific and background signals. These findings establish that sci-mtChIL-seq achieves performance nearly equivalent to that of single-target sci-ChIL-seq, enabling dual-factor localization analysis at single-cell resolution.

Data integration of sci-mtChIL-seq data in myogenesis
Next, using sci-mtChIL-seq, we sought to identify the characteristic epigenomic states associated with each transcriptional state change (Fig. 2a). To do this, we projected multiple sci-mtChIL-seq datasets into a shared space by Harmony22. We employed a MyoD-dependent skeletal muscle differentiation model23. In this system, NIH3T3 cells harboring TetOn3G/MyoD were induced to express MyoD by adding doxycycline (Dox) (Fig. 2b, Supplementary Fig. 2a). Immunostaining revealed that over 70% of cells became MyoD-positive 24 h after the addition of doxycycline [Dox(+): D0h], while the differentiation marker Myog was detected in ~17% of cells at 24 h (D24h) and ~60% at 72 h (D72h) following induction in low-serum medium (Supplementary Fig. 2b), indicating differentiation efficiency of approximately 60%. Here, sci-mtChIL-seq was performed at four time points (before MyoD expression, D0h, D24h, and D72h). Epigenomic factors analyzed simultaneously with RNAPII were Brg1, H3K4me3, H3K27ac, and H3K27me3.

To evaluate whether RNAPII data captured changes of cell state during differentiation, we clustered and visualized single-cell RNAPII counts (Fig. 2c). Cells exhibited a temporal distribution from MyoD(−) to D72h, consistent with progressive differentiation. Regarding the character of each cluster, the composition of the cells in each cluster showed an increasing proportion of cells from the later stage in the order of clusters 1, 2, 4 and 3, suggesting that each cluster reflects a particular stage of differentiation (Fig. 2d). Furthermore, these clusters included cells derived from different antibody-pair datasets without apparent segregation, indicating successful removal of batch effects (Fig. 2d). To characterize these clusters transcriptionally, we identified differentially expressed genes (DEGs) based on RNAPII occupancy. For each cluster, the top 50 DEGs were extracted in order of increasing variance, and the level of RNAPII accumulation at each gene was visualized in the 200 cells with the highest total counts (Fig. 2e). Clusters 1, 2, and 4 showed common gene transcriptional activity in some DEGs, while cluster 3 showed a similar expression pattern only with cluster 4. UMAP visualization of representative DEGs showed cluster-specific expression: cell cycle genes (Top2a, Mcm6) in cluster 1; fibroblast markers (Ptn, Fn1) in cluster 2; and skeletal muscle genes (Cdh15, Ncam1, Myh3, Mef2a) in clusters 3 and 4 (Fig. 2f). GO enrichment analysis was performed using all detected DEGs (Fig. 2g). GO terms related to the cell cycle and cell proliferation were detected in cluster 1, while those related to skeletal myogenesis were detected in clusters 2 and 4. According to the results of Fig. 2c, d, cluster 1 is an undifferentiated and proliferative cell population, cluster 3 is a myoblast-like differentiated cell population, and cluster 4 is one in the process of such differentiation. Meanwhile, cluster 2 expressed fibroblast-related genes, as well as myogenic genes at low levels, indicating a “reprogramming failure”24. Thus, we were able to identify the cell populations that are expected to emerge in this system.

To investigate the chromatin state of these cell populations, we first counted unique fragments for each gene in each paired target and constructed Gene × Cell count matrices. The regions to be counted were set for each paired target (see Methods). For each cluster, the mean accumulation levels were calculated, and the accumulation levels of the genes used in Fig. 2e were visualized (Fig. 2h). The results revealed that epigenetic factors involved in activation, such as Brg1, H3K4me3, and H3K27ac, showed increased accumulation in the DEGs. Meanwhile, the inhibitory H3K27me3 showed a signal distribution that contrasted with that of activation markers, especially in cluster 1. In clusters 3 and 4, H3K27me3 partially accumulated at DEGs, suggesting that the repressing mechanisms of the genes varied. Next, we visualized the accumulation levels of each epigenomic factor on UMAP representation (Fig. 2i). Cluster-specific accumulation was observed for the activation marker, while H3K27me3 showed relatively heterogeneous accumulation for the entire cluster.

These findings indicate that RNAPII-based data integration via sci-mtChIL-seq reveals coordinated changes in cellular and epigenomic states during muscle differentiation.

Application of sci-mtChIL-seq to embryonic myogenesis
To confirm the applicability of sci-mtChIL-seq to in vivo tissues, we performed epigenomic profiling of mouse embryos in the early stages of myogenesis (Supplementary Fig. 3a). We collected trunk tissues from E10.5 and E11.5 embryos, and limb tissues from E12.5, E13.5, and E14.5 embryos25. For these samples, sci-mtChIL-seq was conducted using antibody sets targeting RNAPII-H3K4me3 and RNAPII-H3K27me3. Using this dataset, we attempted to annotate cell types according to the detected clusters and capture the cell-type-specific epigenomic states.

First, we performed clustering and UMAP visualization using RNAPII counts of the dataset, and the cell types were identified referring to known marker genes (Supplementary Fig. 3b)26,27. The identified cell populations included undifferentiated mesoderm, skeletal muscle, vascular, hematopoietic, and dermal lineages, consistent with the expected cellular composition of developing limb tissues from E10.5 to E14.527. The temporal distribution of these cell types was also analyzed (Supplementary Fig. 3c). Undifferentiated mesoderm, erythrocytes, and hepatocytes, typically absent from the limb27, were predominantly derived from E10.5 and E11.5 embryos. By contrast, limb mesenchyme and connective tissue progenitors were primarily derived from E12.5 onwards28. Next, we examined the mean expression levels of marker genes across clusters (Supplementary Fig. 3d). Although most cell types showed specific expression of each marker gene, the limb mesenchyme progenitor and the connective tissue progenitor showed very similar expression patterns26. Furthermore, by aggregating the fragments for each cell type and visualizing the signals on the marker genes using IGV, RNAPII and H3K4me3 signals showed expected patterns, with RNAPII enriched across gene bodies and H3K4me3 near TSSs (Supplementary Fig. 3e). Conversely, H3K27me3 showed distinct, cell-type-specific patterns of gene silencing, such as its accumulation at the Wnt5a locus in hepatocytes.

These results suggest that sci-mtChIL-seq can characterize heterogeneous chromatin states in vivo at high resolution and is applicable to the analysis of developmental epigenomic regulation across complex embryonic tissues.

Transcription factor dynamics in myogenesis at the single-cell level
Next, we attempted to elucidate the dynamics of multiple transcription factors (TFs) during skeletal muscle differentiation by integrating sci-mtChIL-seq data. The procedure for analysis of binding dynamics in each peak using two sets of sci-mtChIL-seq data is schematically illustrated in Fig. 3a. Through data integration, we inferred the binding dynamics of two TFs along a shared pseudotime. As a model, we used the skeletal muscle differentiation system, as in Fig. 2, and focused on the master TF MyoD29 and its downstream TF, Myog.

The detected clusters were similar to those detected in Fig. 2, and no bias due to a paired target was observed (Fig. 3b). To specifically analyze muscle differentiation, we selected clusters 1, 3, and 4 related to muscle differentiation, after PHATE embedding30 with pseudotime was estimated in the inferred space. Density plots of cells validating pseudotime showed no dataset-dependent batch effects and captured variations reflecting sampling time (Fig. 3c). MyoD(−) cells distributed to the early pseudotime phase, while D72h cells distributed across early to late pseudotime stages (Fig. 3c, bottom), confirming that the inferred pseudotime accurately represented the differentiation trajectory. Next, we grouped cells in a differentiation-state-dependent manner to stably estimate TF dynamics. First, to exclude MyoD(−) and equivalent cells, we estimated the threshold separating MyoD(−) and MyoD(+)(D0h, D24h, D72h) samples and removed cells that did not meet this threshold.

The remaining cells were divided into 18 groups (pseudotime group, PG) based on pseudotime (Fig. 3d). The total number of cells tended to be smaller as the PG progressed, but the proportion of D72h cells increased. Thus, PG also yielded a pattern for the progression of differentiation as pseudotime. The total number of MyoD fragments increased slightly with PG (Fig. 3e). Next, MyoD binding patterns during differentiation were estimated for each peak based on PG (Fig. 3f), and the peaks were subsequently classified into four groups: monotonically decreasing, monotonically increasing, late increasing and transiently increasing (Fig. 3g). GO enrichment analysis of genes overlapping with peaks in the promoter showed that clusters 1 and 2 were predominantly enriched for genes involved in muscle differentiation, whereas cluster 3 was enriched for genes related to transcription (Fig. 3h). These results suggest that MyoD plays a role in inducing myogenic programs up to intermediate stages, while it may continue to regulate the transcriptional state at later stages.

We analyzed RNAPII-Myog in the same manner as used in the RNAPII-MyoD dataset analysis. The results revealed that Myog showed a marked increase in the number of fragments at later stages (Fig. 3i), compared with MyoD, and the majority of peaks showed a monotonically increasing trend in the dynamics of each peak (Fig. 3j, k). In GO enrichment analysis, GO terms related to skeletal muscle formation were observed in peak cluster 1, showing a monotonous increase (Fig. 3l).

Taken together, these results indicate that muscle-related genes are regulated in a sequential manner: initially by MyoD (peak clusters 1 and 4), followed by Myog (peak cluster 1). This supports a model in which MyoD binding precedes Myog recruitment, facilitating the transcription of muscle-specific genes during differentiation31,32. Thus, sci-mtChIL-seq effectively enabled inference of TF binding dynamics, revealing stage-specific regulatory programs underpinning muscle differentiation.

Combinatorial analysis using integrated sci-mtChIL-seq data
To investigate how multiple epigenomic factors, including the analyzed TFs, cooperatively regulate transcription at individual loci, we developed a framework to infer the regulatory order of factor binding events. We describe this strategy as the Epigenome trajectory.

The workflow of the analysis is shown in Fig. 4a and Supplementary Fig. 4a. RNAPII counts from multiple sci-mtChIL-seq datasets—each paired with different epigenomic factors—were embedded into a common low-dimensional space. To impute information from different modalities, we aggregated cells with similar transcriptional states into a meta-cell across the datasets. For each meta-cell, the mean signal for each factor in the binned genomic region was calculated to obtain a (Target × Bin) × Meta-cell matrix. To focus on the regulation of a gene, only the bins overlapping the gene locus of interest (LOI) were extracted. Here, we extended GeneTrajectory33, which estimates the order of gene expression during a biological process using a Gene × Cell matrix from single-cell RNA-seq data. GeneTrajectory was applied to the (Target × LOI) × Meta-cell matrix to estimate the order of accumulation of each target at LOI; Epigenome trajectory. It should be noted that this method only estimates the order of expression, not the time interval or whether there is an increase or decrease in the level of accumulation. With this framework, we analyzed the five sci-mtChIL-seq datasets targeting Brg1, H3K4me3, H3K27ac, MyoD, and Myog simultaneously with RNAPII, to estimate epigenome trajectories (Supplementary Fig. 4b–d).

We first analyzed epigenomic regulation, focusing on the skeletal muscle cell markers Mef2a, Lmod3, and Xirp1, and the non-skeletal muscle cell (cluster 2) marker Ptn (Fig. 4b). For myogenic genes, MyoD involvement commonly appears early in the predicted Epigenome trajectories and predates the appearance of RNAPII and most of the other events. Notably, MyoD binding occurred at different time points for different genes (early for Lmod3, intermediate for Mef2a, intermediate–late for Xirp1), but the sequence of subsequent events was similar between genes. Ptn expressed in cluster 2, the non-muscle cell population, did not show binding of MyoD or Myog, but showed a similar accumulation pattern for the other factors. This suggests that a MyoD-independent regulatory cascade controls the transcription, which accelerates the emergence of cluster 2.

Next, we performed a detailed analysis of Mef2a, a skeletal muscle-specific gene. For Mef2a, each feature was assigned an order along the Epigenome trajectory estimated in Fig. 4b (order), and those with orders 1–5, 1–30, 1–70, 1–90, and 1–126 are shown (Fig. 4c). First, there was accumulation of Brg1 and H3K27ac in the enhancer region on the gene body (order 1–5), followed by MyoD binding in the promoter (order 1–20), accumulation of H3K4me3 (order 1–50), Myog binding (order 1–80) and finally RNAPII elongation (order 1–126). The initial accumulation of Brg1 and H3K27ac occurred in a region overlapping with the distal enhancer defined by ENCODE2,34 (Supplementary Fig. 4e). Furthermore, bulk ATAC-seq results indicated that the same region is accessible at the MyoD(−) stage (Supplementary Fig. 4f). These results suggest that the accessibility of the Mef2a enhancer contributes to transcriptional priming in NIH3T3 cells.

Taken together, these findings support a model in which MyoD binds at enhancers and promoters to promote the formation of active chromatin structures, followed by Myog binding and RNAPII recruitment35 (Fig. 4d). Meanwhile, even though the timing of MyoD binding differs, the downstream regulatory sequence is largely the same, indicating a coordinated epigenomic regulatory mechanism for transcriptional activation.

Methods
Cells
NIH3T3 Tet-On 3G/pCW-MyoD-Blast cells were cultured in DMEM supplemented with 10% Tet-Approved FBS (Takara; 631101). Cells were seeded in 10-cm dishes at 40%–50% confluency (500 K cells per dish) and cultured overnight. The next day, Doxycycline (Nacalai Tesque Inc., 19088-71) was added at a final concentration of 1 μg/ml, and the cells were cultured overnight. Differentiation was induced by changing the culture medium to DMEM supplemented with 2% horse serum (Equitech-Bio, SE30), and cells were collected after 24 and 72 h.

Collection of mouse embryo tissue and cell isolation
Pregnant C57BL/6NCr mice were purchased from Japan SLC, Inc. (Shizuoka, Japan) and maintained under the vendor’s specific pathogen-free conditions, with a 12-h light/dark cycle, ambient temperature of 23–25 °C, and 45–65% relative humidity until embryo dissection. The trunks at E10.5 and E11.5 and the limbs at E12.5, E13.5, and E14.5 were dissected from the embryos and cut into small pieces approximately 1 mm square. For each dataset (for example, E10.5 RNAPII-H3K4me3), tissues were collected from 18–39 embryos. Tissue pieces were suspended in Gey’s Balanced Salt Solution (Sigma; G9779) and transferred to 2 ml low-absorption tubes. Dispase® II (neutral protease, grade II) (Merck; 04942078001) was added to a final concentration of 1.8 U/ml, and the mixture was gently shaken for 10 min at room temperature. After the reaction, the tubes were centrifuged at 1500 × g for 1 min at 4 °C, and the supernatant was discarded. TrypLE™ Express Enzyme (1×), no phenol red (Thermo Fisher Scientific; 12604013), was then added to the pellet, which was gently shaken at 37 °C for 7 min. DNase I (Takara, 2270A) was added at a concentration of 0.014 U/μl, and the cells were gently suspended, washed once with DMEM (Thermo Fisher Scientific, C11885500BT) containing 10% FBS (NICHIREI, 20L00J), resuspended, and passed through a 40 μm cell strainer, after which the cells were counted for sci-mtChIL-seq.

Antibodies
The following antibodies were used: mouse anti-RNAPII S2ph (CMA602, RRID: AB_2819246, clone pc26), rabbit anti-MyoD (Cell Signaling Technology; cat. no. 13812S, RRID: AB_2798320, clone D8G3), rabbit anti-Brg1 (Abcam; cat. no. ab110641, RRID: AB_10861578, clone EPNCIR111A), rabbit anti-Myog (Abcam; cat. no. ab124800, RRID: AB_10971849, clone EPR4789), rabbit anti-H3K4me3 (Abcam; cat. no. ab213224, RRID: AB_2923013, clone EPR20551-225), rabbit anti-H3K27me3 (Merck; cat. no. 07-449, RRID: AB_310624), and rabbit anti-H3K27ac (Active Motif; cat. no. 39133, RRID: AB_2561016). For sci-ChIL probe creation, the following antibodies were used: AffiniPure™ Rat Anti-Mouse IgG (H+L) (Jackson ImmunoResearch Laboratories; cat. no. 415-005-166) and AffiniPure™ Goat Anti-Rabbit IgG (H+L) (Jackson ImmunoResearch Laboratories; cat. no. 111-005-144).

Immunocytochemistry
Cells were plated on ibidi 96-well plates, washed twice with phosphate-buffered saline (PBS), fixed with 1% paraformaldehyde (PFA, Electron Microscopy Sciences, 15710) in PBS, permeabilized with 1% Triton X-100 in PBS, and washed twice with PBS. After incubation with Blocking One-P solution (Nacalai Tesque Inc.) for 15 min, the cells were incubated with rabbit anti-Myog (1:400) or rat anti-MyoD (Kishida; cat. no. 380-00006, RRID: AB_2793421, clone 5F11, 1:200) diluted with 0.1× Blocking One-P/PBS (mix of two antibodies for sci-mtChIL-seq) Blocking One-P in PBS for 2 h at room temperature. The cells were then washed three times with PBS, supplemented with Alexa546-conjugated goat anti-rabbit or anti-rat antibodies (1:2000; Thermo Fisher Scientific) and bisbenzimide H33342 fluorochrome trihydrochloride (Hoechst) (1:2000; Nacalai Tesque Inc.), followed by incubation for 30 min at room temperature, and then the cells were washed three more times with PBS. Images were obtained with a fluorescence microscope (BZ-X710; Keyence). Co-localization was evaluated using the CellInsight CX5 HCS Platform (Thermo Fisher Scientific, CX51110) from the nine different field images.

Bulk RNA barcoding and sequencing (BRB-seq) for purified RNA
Cells were cultured in a 6-well plate (Cat. No. TR5000 TrueLine Cell Culture Plate). Doxycycline (Dox) was added to the culture medium at a final concentration of 1 μg/ml for 24 h. To induce differentiation, the medium was replaced with DMEM supplemented with 2% horse serum (Equitech-Bio, SE30), and cells were harvested at 24 and 72 h post-induction. Total RNA was extracted using TRIzol™ Reagent (Ambion) and purified using the Direct-zol RNA Microprep Kit (ZYMO RESEARCH). RNA quality was assessed using the TapeStation 4200 system (Agilent Technologies), and concentration was measured using the Qubit™ RNA High Sensitivity Assay Kit (Invitrogen). BRB-seq was performed using 40 ng of total RNA per sample. The protocol followed BRB-seq v2.039, with one modification in cDNA synthesis: reverse transcription was performed using CEL-Seq2 RT primers40 instead of the original BRB-seq primers. The quality of the cDNA was assessed using TapeStation and quantitative PCR. Sequencing was performed on an Illumina NovaSeq 6000 system using the NovaSeq 6000 SP Reagent Kit v1.5 (100 cycles) using 19 cycles for Read 1 and 51 cycles for Read 2.

ATAC-seq
Cells were grown in 15 cm dishes. Doxycycline was added at a final concentration of 1 μg/ml for 24 h. Differentiation was induced by changing the medium to DMEM supplemented with 2% horse serum (Equitech-Bio, SE30), and cells were harvested after 24 and 72 h. Approximately 50,000 cells were washed with cold PBS and lysed with cold lysis buffer (1× NP40, 1× Tween-20, 1× Digitonin). After incubation on ice for 3 min, the cells were centrifuged at 500 × g for 10 min at 4 °C. The resulting pellet was resuspended in the transposition reaction mix containing 1× TAPS+DMF buffer, 100 nM Tn5 transposase, 0.35× PBS, 1× digitonin, and 1× Tween-20 in nuclease-free water. The reaction was incubated at 37 °C for 30 min and immediately placed on ice. DNA was purified using the ZYMO DNA Clean & Concentrator-5 Kit (D4014, Zymo Research). The resulting library fragments were PCR amplified using custom primers for 5–6 cycles and purified again using the same kit. The libraries were sequenced on an Illumina NovaSeq 6000 system using the NovaSeq 6000 SP Reagent Kit v1.5 (100 cycles), with 51 cycles each for Read 1 and Read 2.

Production of sci-ChIL probe
The sci-ChIL probe was prepared based on our previously established protocol13,14. Oligonucleotides were synthesized by Eurofins Genomics and purified by high-performance liquid chromatography. Rat anti-mouse and Goat anti-rabbit IgG antibodies (Jackson ImmunoResearch; cat. no. 415-005-166 and 111-005-144) were conjugated with DBCO-PEG5-NHS ester (Click Chemistry Tools) using a Cu2+-free strain-promoted alkyne–azide cycloaddition reaction. The reaction mixture contained 100 μg antibody in 100 mM NaHCO3 (pH 8.3) with a 10-fold molar excess of DBCO-PEG5-NHS ester, incubated for 1 h at room temperature. The conjugate was desalted using a PD MiniTrap G-25 column (Cytiva), concentrated with an Amicon Ultra-0.5 NMWL 10 kDa centrifugal filter (Merck Millipore), and adjusted to 1 μg μl−1 in PBS. Double-stranded ChIL DNA was prepared by annealing two complementary oligonucleotides (95 °C for 5 min, then cooled at −0.1 °C s−1 to 20 °C) in 10 mM Tris-HCl (pH 7.4), 1 mM EDTA, and 100 mM NaCl. Antibody–DNA conjugation was performed by mixing 75 μl of DBCO-antibody (1 μg μl−1) with 100 μl of 10 μM DNA and incubating for 1 week at 4 °C. Unreacted DNA was removed by five rounds of filtration through a Microcon Ultracel DNA Fast Flow Membrane (Merck Millipore). All oligonucleotide sequences used, including those for cell barcoding and antibody conjugation, are provided in Supplementary Data 1.

sci-ChIL-seq and sci-mtChIL-seq protocol
sci-ChIL-seq and sci-mtChIL-seq were performed the developed based on our previously described methods12–14. Each process is described below.

Cell surface biotinylating. To 1×107 cell pellets, 1 ml of ice-cold 1 mg/ml Biotin-SS-Sulfo-OSu (DOJINDO) in PBS to suspend the pellets, followed by rotation at 4 °C for 30 min, washing twice with quench buffer [50 mM Tris-HCl (pH 8.0), 0.1 mM EDTA, 150 mM NaCl], and washing twice with PBS. Next, 1% PFA/PBS was added to the cell surface-biotinylated cells, the cells were fixed at room temperature for 5 min, and 1 M glycine was added, followed by incubation at room temperature for 5 min. After two washes with PBS, the cells were suspended in PBS/0.1% BSA/0.05% Tween-20 and dispensed into 96-well Pierce™ NeutrAvidin™-Coated High-Capacity Plates, Clear, 8-Well Strip (Thermo Fisher Scientific), at 50 μl/well (≤10,000 cells). The wells were left overnight at 4 °C with slow shaking and then washed twice with PBS.

Immunostaining. Cells were permeabilized with 1% Triton X-100 in PBS for 20 min and blocked with Blocking One-P (Nacalai Tesque Inc.) for 20 min at room temperature. Then, 2 μg/ml primary antibody solution in 0.1× Blocking One-P/PBS (mix of two antibodies for sci-mtChIL-seq) was added to the wells, which were shaken slowly at 4 °C for 6 h. After three washes (5 min each at 4 °C) with PBS, 2 μg/ml sci-ChIL probe solution was added to 0.1× Blocking One-P and 0.5 M NaCl/PBS (for sci-mtChIL-seq, the sci-ChIL probe combination is shown in Supplementary Data 2), and the mixture was shaken slowly overnight at 4 °C.

ChIL reaction. Tn5 (90 ng of purified Tn5) in 1× Tn5 dialysis buffer (50 mM HEPES-KOH, pH 7.2, 0.1 M NaCl, 0.1 mM EDTA, 1 mM dithiothreitol, 0.1% Triton X-100, and 10% glycerol) was added to each well and incubated for 10 min at room temperature. The Tn5 transposase was produced in-house following a previously described protocol41. MEDS-B (10 μM, 0.1 μl/well) in 1× Tn5 dialysis buffer was then added and allowed to react for 1 h at room temperature. After three washes with PBS (5 min at RT), 1× TAPS-DMF buffer (10 mM TAPS-NaOH, pH 8.5, 5 mM MgCl2, and 10% N,N-dimethylformamide) was added and reacted at 37 °C for 1 h for the Tn5-mediated integration reaction. The reaction was then stopped with 50 mM EDTA at 50 °C for 30 min, the mixture was washed three times with PBS, and then a fill-in solution (200 U per well T4 DNA ligase, NEB, 1.5 U T4 DNA Polymerase, NEB, in 1× T4 DNA Ligase Reaction buffer, NEB) was added, followed by incubation at room temperature for 30 min. The reaction was stopped with 50 mM EDTA at 50 °C for 30 min, followed by washing three times with PBS. After the addition of Hoechst to the wells and nuclear staining of the cells, Accumax (Nacalai Tesque Inc., 17087-54) for SP2 cells or a solution of Accumax: Accutase (Nacalai Tesque Inc., 12679-54) = 1:3 for embryonic cells and NIH3T3 cells was added and allowed to settle for 10 min at room temperature. Additionally, 1 M DTT (Nacalai Tesque Inc.) was added, and allowed to react at 37 °C for 30 min to allow the cells to detach from the plate. The entire supernatant was collected in a 15 ml tube and centrifuged to precipitate the cells. Approximately half of the volume of the supernatant was removed, a final concentration of 0.1% Tween-20/2% BSA in PBS was added, and the cells were resuspended. Sorting was performed on a SONY SH800. Cell solution was filtered through a 40 μm filter before sorting, and Hoechst-positive cells were sorted to 25 cells/well in a 96-well plate containing IVT solution (2 U/μl per well; TOYOBO). IVT reactions were then incubated at 37 °C for 3 h. RNA samples were treated with DNase I (2.5 U per well; TaKaRa) for 30 min at 37 °C. The treated RNA samples were purified by adding 1.8× volume of Agencourt RNAClean XP beads (Beckman Coulter). They were then incubated with the beads for 5 min at room temperature, washed twice gently in 70% ethanol, and eluted in 7 μl of nuclease-free-water. Reverse transcription was carried out in accordance with the protocol of the SMART-Seq v4 Ultra Low Input RNA Kit (Clontech/Takara) with minor modifications, but all reagents were prepared separately. In each well of a 96-well plate, 7 μl of purified RNA were added in 1.4 μl of reaction buffer [0.655 μl of 10× Lysis Buffer, 0.7 μl of Read 2 primer (12 μM), 0.035 μl of RNase inhibitor]. Samples were incubated at 72 °C in a preheated, hot-lid thermal cycler for 3 min. Then, the samples were supplemented with 2.8 μl of 5× First-Strand Buffer, 1 μl of Advantage Ultra Pure PCR Deoxynucleotide Mix, 1 μl of DTT, 0.25 μl of RNase Inhibitor, and 1 μl of SMART Scribe Reverse Transcriptase. The reverse transcription was performed in a thermal cycler with the following program: Step 1: 42 °C × 90 min, followed by Step 2: 70 °C for 30 min and holding at 4 °C. After the reaction, the samples were used for final PCR with the following reaction mixture: 14.45 μl of template cDNA, 15 μl of 2× SeqAmp PCR Buffer, 1 μl of unique index primer (10 μM, Supplementary Data 3), and 0.5 μl of SeqAmp DNA Polymerase. PCR was performed at 95 °C for 1 min; 16–18 cycles of 98 °C for 10 s, 65 °C for 30 s, and 68 °C for 3 min; and final extension at 72 °C for 10 min. RNase A (5 μg per well; Nacalai Tesque Inc.) was then added for 30 min at 37 °C to remove any remaining RNA. The 96-well amplified DNA was combined in a 15 ml tube and purified, and then size-selected using 0.55× and 0.9× two-sided purification. The sample was supplemented with 0.55× volume of Axygen AxyPrep MAG PCR Clean-Up (CORNING), followed by incubation with beads for 5 min at room temperature, after which the supernatant was collected. Then, the supernatant was supplemented with 0.35× volume of beads for 5 min at room temperature, and the supernatant was discarded. The beads were washed in 80% EtOH and eluted in 100 μl of 10 mM Tris, pH 8.0. The size selection was then repeated with the same percentage of beads. This process removed all sequences except those of approximately 200–800 bp in size. Libraries were quality-checked and quantified by Tapestation (Agilent Technologies) and qPCR (Bio-Rad). The libraries were sequenced on NovaSeq 6000 using custom primers as follows: i5: 35 bp, i7: 8 bp, read 1: 46 bp, and read 2: 46 bp.

MEDS-B
Tn5MErev: CTGTCTCTTATACACATCT
Tn5ME-B: GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG
Read 2 primer: GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGNNNNNNNNNNNN
Custom primers
Read 1 Seq Primer (Extended ME-A): GCGATCGAGGACGGCAGATGTGTATAAGAGACAG
Custom index 2 primer: CTGTCTCTTATACACATCTGCCGTCCTCGATCGC

Single-cell data analysis
Read mapping and processing. To obtain the sequences of both the reads and dual indices, run the following command: “bcl2fastq --create-fastq-for-index -R bclfile_directory.” Cell barcodes were made by modifying the index fastq files by seqkit (version 0.16.1)42. Cell barcodes were added to the read names using umi_tools (version 1.1.1)43. Adapter sequences were trimmed using TrimGalore (version 0.6.10) with the default options. Paired reads were mapped to mm10 by Bowtie2 (version 2.3.5.1)44 with the default options. BAM files were created, sorted, and indexed by SAMtools (version 1.18)45. Paired reads mapped to autosomes and sex chromosomes were obtained. To obtain unique reads, PCR duplicates that had the same cell barcode and were mapped to the same site were excluded. Unique reads with the same cell barcode and with the same start position are expected to be caused by random primers in IVT. Unique fragments were obtained by removing the reads. The number of obtained single-cell data and the average read counts per cell were provided in Supplementary Data 4.

Quality control. Out of the expected 96×96 barcodes, up to 2400 barcodes were assigned to a cell in the experimental design in this study. Therefore, barcodes with a low number of fragments, which were not assigned to any cell, and detected by chance. To exclude these barcodes, thresholds for the number of fragments were set for each dataset. To calculate appropriate thresholds, the distributions of the number of fragments for assigned and unassigned barcodes were estimated, respectively, using the EM algorithm implemented in the R package mixtools (v2.0.0)46. Thresholds were set when the cumulative probability of the assigned barcode distribution reached 0.01.

Additionally, to evaluate the risk of doublets (collisions of the cell barcodes) in the combinatorial indexing procedure, we created the n×m contingency table of the number of reads obtained X=(xij) for all combinations of primary (T5) and secondary index (i5+i7). We assumed that (a) the first stage of the T5 indexing represents a non-uniform distribution with the barcode assignment probabilities: pi for i=1, 2, …, n, while (b) the probabilities are uniform in the second stage of i5+i7 indexing: qj = 1/m for j=1, 2, …, m. The assumptions (a) and (b) give the expected number of cells with a combination of the barcodes (i,j) as λij = Cpi/m, where C is the total number of cells used. Since the mean parameter λ is constant λij = λi over j (during the second indexing), m observations in each ith row of X can be used as m i.i.d. samples from a Poisson distribution with parameter λi. Note that we can only observe that the number of assigned cells k in is zero (k=0; no cells assigned to the barcode) or non-zero (k>0; singlet, doublet, or more). The probability of k=0 is f(k=0; λi) = e−λi in the Poisson distribution; therefore, λi was estimated by λ^i = −log zi/m where zi is the number of unobserved barcodes in the ith row of X. The zi can be the number of barcodes for which the number of reads xij was less than a threshold in the ith row. Given a barcode (i, *), the probability r = f(k≥2; λ^i) can be used as a risk index of doublets. We excluded cells with r > 0.2 and a unique read count greater than twice the median.

To correct Tn5 bias, the regions detected as peaks by MACS247 for all datasets were blacklisted, and fragments overlapping the regions were excluded.

Quantification. Unique fragments were counted for each gene. For the reference genes, protein-coding genes and their loci were obtained from EnsDb.Mmusculus.v79. For genes with variants, the longest one is used as a representative.

The region for counting fragments was set as follows: RNAPII and H3K27me3; from TSS+1 kbp to TES+1 kbp, H3K4me3 and Brg1; TSS±5 kbp, H3K27ac; gene body.

Dimensional reduction. Dimensional reduction and clustering were performed following Seurat (version 5.3.0)48. Batch effects were removed by Harmony22 using principal components (PCs). PCs for downstream analysis were set based on an elbow plot. For all datasets, the correlation coefficients between Harmony components 1 and read count per cell were nearly |1|, so the components were excluded for the following analysis. Clustering was performed with the Seurat functions FindNeighbors and FindClusters. In FindClusters(), a Leiden algorithm was used.

Cell type annotation. Based on marker genes used by He et al. 27 and Cao et al. 26 for cell type annotation in scRNA-seq datasets, the clusters in sci-mtChIL-seq were evaluated using scType49. The cell types of the clusters were identified with reference to scType scores.

Integrative Genomics Viewer (IGV). To visualize the aggregated signals distribution, the BigWig files were created by deeptools bamCoverage (version 3.5.4)50 for fragments of each dataset with the option --normalizeUsing CPM, –centerReads, -binSize 100, and --smoothLength 1000. For H3K27me3 fragments, the following values were used: -binSize 1000 --smoothLength 10000.

Gene Ontology (GO) enrichment analysis. GO enrichment analysis for biological processing was performed with the R package clusterProfiler (version 4.6.2)51 and org.Mm.eg.db.

Pseudotime analysis. The trajectory was inferred using the R package phateR (version 1.0.7)30 within the same principal component analysis space used for UMAP. Cells not involved in the biological process of interest, specifically non-muscle lineage cells, were excluded from the trajectory inference. Pseudotime was estimated using slingshot (version 2.6.0)52, based on the PHATE space. MyoD-negative cells were removed using a threshold estimated to distinguish between MyoD-positive and MyoD-negative samples.

The pseudotime was divided into 20 equally sized segments, referred to as PGs. PGs containing a small number of cells were merged with neighboring groups. For each PG, cells were aggregated, and fragment counts were obtained for each TF peak from ENCODE (MyoD: ENCFF814ERT, ENCFF138LQJ, ENCFF779ENE, ENCFF991ZUO; Myog: ENCFF277SLZ, ENCFF695DBU, ENCFF659JQU, ENCFF766TEE). These counts were used as estimates of TF binding levels.

To infer binding dynamics at each peak, binding levels were smoothed across PGs using a generalized additive model53, and the predicted values were scaled to the range [0, 1]. Peak clusters were identified via k-means clustering based on their dynamic profiles, with the number of clusters determined using silhouette analysis. Genes associated with promoter-proximal peaks (from −2 kbp to the TSS) were used for GO enrichment analysis.

Epigenome trajectory. Epigenomic regulation was inferred using the GeneTrajectory method33. Originally, this method estimates the sequential order of gene expression within a specific biological process, based on a low-dimensional representation of the Wasserstein distance between genes. The Wasserstein distance is computed pairwise between distributions of gene expression levels on a cell–cell graph, with the transportation cost defined by the shortest path. To reduce computational cost, a coarse-grained version of the cell–cell graph was used for the Wasserstein distance calculation.

In this study, we modified the processes of cell–cell graph construction and Wasserstein distance calculation to accommodate the sci-mtChIL-seq datasets. The cell–cell graph was constructed using the latent representation of RNAPII counts for binned coding regions (bin width=5 kbp), estimated by PeakVI54 to correct for batch effects and to estimate denoised counts. The count matrices for modalities other than RNAPII were also processed using PeakVI. A total of 1000 meta-cells were defined from the cell–cell graph using k-means clustering. For the coarse-graining step, denoised counts for all modalities were aggregated at the meta-cell level. Wasserstein distances were calculated among regions within the gene LOI to infer gene-specific epigenomic regulation.

For visualization of the results, a principal curve was inferred, and features were projected onto this curve55. The distributions of regulatory factors along the projected coordinates were estimated using the R package ggridges. The expression order was also determined based on these coordinates.

Bulk data analysis
Bulk RNA-seq (BRB-seq) data analysis. Read 1 (barcode read) was extracted by using UMI-tools with the following command “umi_tools extract -I read1.fastq --read2-in=read2.fastq --bc-pattern=NNNNNNNNNNCCCCCCCCC --read2-stdout ”. Adapter sequences and low-quality sequences were removed, and read lengths below 20 bp were discarded by using TrimGalore, and reads were mapped to the GRCm38 reference using HISAT2 (2.2.1)56. Read counts for each gene and each sample were obtained by summarizing the BAM file, which was produced by featureCounts (2.0.1)57.

Bulk ATAC-seq data analysis. Paired reads were mapped to the mm10 reference using Bowtie2. Alignments that were unpaired, discordant, or had fragment lengths exceeding 1000 bp were discarded. BAM files were created, sorted, and indexed by SAMtools. Bigwig files were created by the following command “bamCoverage --binSize 100 --centerReads --normalizeUsing CPM -b BAM.bam -o 100PE.CPM.bw”.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This study by Harada et al. investigates the regulatory mechanisms underlying myogenic differentiation during mouse embryonic development using single-cell multi-omics techniques. The authors profiled chromatin interactions and transcriptional activities in myogenic progenitor cells across key developmental stages. They employed a sci-mtChIL-seq methodologies to map chromatin accessibility, histone modifications, and transcription factor binding dynamics with high spatial and temporal resolution. Their findings reveal dynamic changes in chromatin architecture and transcriptional regulation orchestrated by key regulatory factors such as MyoD and Myog. Moreover, the study highlights the utility of single-cell multi-omics approaches in unraveling complex gene regulatory networks during embryonic development.\nSpecific comments are as follows:\n1. All conclusions drawn in the manuscript were not validated experimentally, which significantly undermines the strength of their claims.\n2. The method description was not informative. This is critical because the conclusion could lead to overfitting.\n3. Page 4 line 4 It is not clear how to set cut-off values and cell-derived indices.\n4. Page 4 line 8 “The genomic localization patterns of RNAPII from RNAPII-H3K4me3 and RNAPII-H3K27me3 were roughly similar between stages.” How can this result be explained, and does it indicate that the localization of PolII is not related to the levels of H3K4me3 and H3K27me3?\n5. Page 5 line 4 What is stage-dependent bias?\n6. Fig 3e-f MyoD are inconsistent between Lineage1 and Lineage2.\nOther comments:\nThere are numerous sections of the text that require major editing for clarity, accuracy, and to provide enough information.\nBelow are a few examples:\nFig 1i left RNAPIII\nFig 4f RNAII\nExtended Data Fig 3a left tH3K4me3", "reviewer_2": "Harada and coauthors present multi-target ChIL-seq (mtChIL-seq) to simultaneously analyze the binding of TFs, RNAPII, and other molecules to chromatins and the high-throughput extension of this chemistry to develop single-cell combinatorial indexing mtChIL-seq (sci-mtChIL-seq) that can indirectly trace cell lineages with RNAPII and evaluate epigenomic status. They performed sci-mtChIL-seq on the lineage of embryonic myogenesis of mouse limbs to classify cell types based on the binding state of RNAPII and simultaneously profile changes in the epigenomic state. They also analyzed the role of TFs in cell lineage determination by invetigating MyoD-dependent epigenomic state development using in vitro models. Overall, this manuscripts presents some quite interesting data. The development of sci-mtChIL-seq itself in this work is largely based on the ChIL–seq method developed by the same group, and provides a single-cell ChIL-seq profiling solution with the split-pool barcoding strategy. There are some technical comments.\nThe authors provide solutions to both single-target sci-ChIL-seq and the multi-target sci-mtChIL-seq. But the first section of the manuscript only covers the single-target sci-ChIL-seq. Neither the context or the figures cover any details regarding how multiple targets can be profiled at the same time. Did the authors incubate primary antibodies that target multiple antigens at the same time? In the multi-target method (for example RNAPII and H3K27me3), are the sequencing reads aggregates of RNAPII and histone signals, or the authors can stratify/demultiplex the RNAPII and H3K27me3 profiles?\nAfter the bulk ChIL–seq was described in 2018, methods like CUT&RUN and CUT&Tag turn to be the more accurate and robust for histone or epigenetic factor binding profiling. Can the authors compare the data throughput, sensitivity, specificity, and performance with these related methods?\nMore details of this protocol should be described, such as why the authors chose to distribute a total of 25 cells per well in split-pool barcoding. Is FACS sorting required for this procedure? It is recommended to mention key procedures such as cell fixation, permeabilization and blocking in the main text to help readers get a better understanding of this technology.\nThe authors mention they use homemade Tn5. Please describe the Tn5 preparation in method and add relevant citations for the protocol.", "reviewer_3": "Harada et al develop a new assay to probe RNA Pol II binding and epigenomic factors in single cells. They then use this tool to evaluate binding dynamics during myogenesis within the developing murine limb. While this is an exciting tool, it’s not clear what new information is included in this manuscript. The results are not clearly articulated and not discussed in the context of what is known.\nExamples of text that are unclear or that require context:\n-The text accommodating Fig. 1i-j says, “..association between cell populations represented by RNAPII binding and histone modifications during limb development.” But what exactly does that tell you? Also, “..the cell population represented by H3K27me3 was..” How do you define large? And why is that significant? These are just two examples of a common theme from the manuscript.\n-Regarding the MyoD binding analyses, what does Fig. 3D show? The text first says a motif enrichment was performed but then it says that the ‘myocyte lineage had the greater deviation score.” What is a deviation score?\n-The authors also reference the term lineage tracing multiple times but as described, their assay does not evaluate lineage as it is traditionally defined. Trajectory inference allows you to infer fate but it does not establish lineage identity."}, "author_response": "We sincerely thank the reviewers for their constructive and insightful feedback. In response, we\nhave undertaken a substantial revision of the manuscript, including the following: experimental\nadditions, major restructuring of the figures, and a complete reorganization of the manuscript to\nenhance clarity and rigor. This revision was driven by our recognition—based on the reviewers’\ninput—that the original submission lacked a compelling scientific focus and that the\ndemonstration of the utility of the developed method was insufficient.\nOver the past year, we have developed a more integrated experimental and analytical framework\nthat better highlights the capabilities of sci-mtChIL-seq. As a result, we have removed parts of\nthe original manuscript that lacked interpretability or reproducibility, and incorporated new data\nand analyses that emphasize a unified, gene-locus-level model of chromatin regulation.\nWhile the revised version may appear substantially different, it represents a conceptual and\ntechnical evolution directly guided by the reviewers’ input. The revised version retains only seven\npanels from the previous version (see Fig. 1a, Supplementary Fig. 1a–d, and Supplementary Fig.\n2a and b), and shows only 32% textual similarity overall. We view the revised version as a direct\nand necessary evolution of the original study. While the initial consideration was to submit this\nas a new manuscript, we felt it was important to honor the reviewer feedback and editorial process\nby presenting the improved version within the same review cycle.\nReviewers' comments:\nReviewer #1 (Remarks to the Author):\nThis study by Harada et al. investigates the regulatory mechanisms underlying myogenic\ndifferentiation during mouse embryonic development using single-cell multi-omics techniques.\nThe authors profiled chromatin interactions and transcriptional activities in myogenic progenitor\ncells across key developmental stages. They employed a sci-mtChIL-seq methodologies to map\nchromatin accessibility, histone modifications, and transcription factor binding dynamics with\nhigh spatial and temporal resolution. Their findings reveal dynamic changes in chromatin\narchitecture and transcriptional regulation orchestrated by key regulatory factors such as MyoD\nand Myog. Moreover, the study highlights the utility of single-cell multi-omics approaches in\nunraveling complex gene regulatory networks during embryonic development.\nSpecific comments are as follows:\n1. All conclusions drawn in the manuscript were not validated experimentally, which significantly\nundermines the strength of their claims.\nWe thank the reviewer for this important point. In the revised manuscript, we addressed this\nconcern by demonstrating that our analytical framework is capable of recapitulating well-\nestablished, experimentally validated regulatory phenomena. Specifically:\n In Fig. 3, our data successfully reproduce the known sequential binding of MyoD followed\nby Myog during myogenic differentiation, consistent with prior perturbation-based studies\n(Cao et al., EMBO J. 2006, PMID: 16437161; Ohkawa et al., EMBO J. 2006, PMID:\n16424906).\n In Fig. 4, our method reconstructs the temporal order of chromatin regulator recruitment at\nmuscle gene loci, in line with established findings on regulatory coordination in muscle genes\n(Hernández-Hernández et al., Front. Genet. 2020, PMID: 33193700; Blum et al., Genes Dev.\n2012, PMID: 23249738).\nThese validations demonstrate that our framework can recover biologically meaningful\nregulatory cascades at single-cell resolution, even in the absence of new perturbation experiments.\nIn the original manuscript, we had focused on demonstrating the utility of our method through\ntechnically consistent patterns derived from known regulatory relationships. However, we now\nrecognize that certain interpretations lacked sufficient biological context or direct support. As a\nresult, we have substantially revised the manuscript, removing under-supported analyses and\nincorporating additional layers of evidence where appropriate.\nFor example, in Supplementary Fig. 5e–f, we show that at the Mef2a locus, the earliest\naccumulation of chromatin modifiers (Brg1 and H3K27ac) occurs at a distal enhancer annotated\nin the ENCODE database, prior to MyoD binding at the promoter, H3K4me3 deposition, Myog\nbinding, and RNAPII elongation. This sequential pattern suggests a temporally ordered cascade\nof regulatory events.\nFurthermore, bulk ATAC-seq data confirm that this distal enhancer region is already accessible\nat the undifferentiated (MyoD–) stage, suggesting transcriptional priming through pre-established\nchromatin accessibility.\nTogether, these findings support a model in which MyoD sequentially engages enhancer and\npromoter regions to promote active chromatin formation, followed by Myog and RNAPII\nrecruitment. The reproducibility of these patterns, their consistency with known enhancer biology,\nand alignment with orthogonal data sources provide strong technical and biological validation of\nour conclusions.\n2. The method description was not informative. This is critical because the conclusion could lead\nto overfitting.\nWe agree with the reviewer that insufficient methodological clarity may weaken the\ninterpretation of our results, particularly with regard to the risk of overfitting. In the revised\nmanuscript, we have clarified the rationale and implementation of trajectory-based visualizations,\nremoved figures that were difficult to interpret due to dataset-specific pseudotime assignments,\nand made the complete analysis pipeline publicly available to ensure transparency and\nreproducibility.\nSpecifically, we inferd that the reviewer’s concern may relate to the trend estimation in the\noriginal Fig. 3e,f and Fig. 4c,g, where we plotted changes in factor accumulation over pseudotime\nusing regression curves. In particular, Fig. 4c,g showed accumulation levels averaged within\npseudotime windows to reduce noise and enhance robustness.\nTo mitigate overfitting, we used LOESS (locally estimated scatterplot smoothing) with a large\nspan parameter, which limits model flexibility. While this reduces temporal resolution, it helps\nprevent the regression from capturing spurious fluctuations. Thus, we believe the risk of\noverfitting in these visualizations was minimal.\nNevertheless, we acknowledge that interpreting trends was complicated by the use of\nindependently computed pseudotime for each dataset, which limited direct comparisons. To\naddress this, we have removed the original panels and, in the revised Fig. 4, present an alternative\nanalysis that integrates multiple datasets using RNAPII as a shared reference. This analysis\nevaluates factor accumulation at individual gene loci, offering a more intuitive and interpretable\nview of regulatory order without relying on pseudotime inference or regression-based smoothing.\nFurthermore, we have made the full analysis pipeline—including data preprocessing,\npseudotime inference, GAM smoothing parameters, and visualization scripts—publicly available\nvia GitHub (see Code availability section). While the revised text outlines the core methodology,\nadditional technical details are available in the repository to ensure full transparency.\nWe thank the reviewer for highlighting this important issue and prompting us to refine the\npresentation and interpretation of our trajectory-based analyses.\n3. Page 4 line 4 It is not clear how to set cut-off values and cell-derived indices.\nWe agree with the reviewer that the method for defining cut-off thresholds and identifying cell-\nderived indices should be clearly explained. In the revised manuscript, we now describe how these\nthresholds were estimated in a data-driven and dataset-specific manner using a statistical model.\nSpecifically, for barcode demultiplexing, we used a combinatorial split-pool indexing strategy,\nin which unique cell identities were assigned via multiple rounds of barcoding. To remove\nbarcodes with low read counts likely derived from ambient background or non-cell-containing\nindices, we used a two-component Gaussian mixture model and estimated the threshold values as\nthe intersection point between the two distributions. Additionally, we inferred doublet probability\nbased on the frequency of occurrence per index and exclude potential doublets or multiplets.\n4. Page 4 line 8 “The genomic localization patterns of RNAPII from RNAPII-H3K4me3 and\nRNAPII-H3K27me3 were roughly similar between stages.” How can this result be explained, and\ndoes it indicate that the localization of PolII is not related to the levels of H3K4me3 and\nH3K27me3?\nWe thank the reviewer for pointing out the ambiguity in our original statement. In response, we\nhave removed the sentence in question from the revised manuscript to avoid misinterpretation.\nIn the revised version, we instead emphasize that RNAPII signals from both RNAPII-H3K4me3\nand RNAPII-H3K27me3 datasets were successfully integrated into a common space using\nHarmony, and showed consistent clustering without major batch-specific segregation\n(Supplementary Fig. 3b). This supports the conclusion that RNAPII localization reflects shared\ntranscriptional activity across stages, while also being compatible with the chromatin context (e.g.,\nactive or poised states).\nWe hope this clarification resolves the concern and makes our interpretation more precise.\n5. Page 5 line 4 What is stage-dependent bias?\nWe thank the reviewer for this important question. In the original manuscript, the term “stage-\ndependent bias” was used to describe an observation in the H3K27me3 dataset, where cell\nclustering patterns corresponded more strongly to developmental stage than to cell type identity.\nHowever, we now recognize that this terminology may have been ambiguous and potentially\nmisleading.\nSpecifically, in the H3K27me3 dataset derived from embryonic limb tissues, UMAP\nvisualization showed the cluster structure that aligned closely with embryonic day (e.g., E10.5 vs.\nE12.5), rather than with transcriptionally distinct cell populations. This pattern likely arises from\nbroad H3K27me3-marked heterochromatin domains, which can be partially shared across\nlineages but shift globally during developmental progression. Thus, the stage-dependent pattern\nreflects a biological phenomenon related to chromatin state transitions, not a technical artifact.\nHowever, we agree with the reviewer that the term “bias” may unintentionally suggest a\nconfounding or technical error. To avoid this confusion, we have removed the term and its\nassociated interpretation from the revised manuscript. This decision was made not due to any data\nquality issue, but to improve clarity and avoid misinterpretation. We emphasize that the embryo\ndatasets remain a core component of our analysis, particularly for RNAPII-based cell state\nclassification and trajectory inference. The removal of this terminology does not reflect a\ndismissal of the embryo data, but rather a refinement of how we present their biological\ninterpretation.\nWe thank the reviewer again for drawing our attention to the potential ambiguity in our wording.\n6. Fig 3e-f MyoD are inconsistent between Lineage1 and Lineage2.\nWe thank the reviewer for this valuable comment. We have addressed this issue by removing the\noriginal MyoD binding analysis in embryonic samples from the revised manuscript.\nAs we noted in the first submission, the number of MyoD-expressing cells captured from\nembryonic tissues was limited, making reliable trajectory reconstruction and quantitative\ncomparison across lineages technically difficult. Moreover, MyoD expression is known to initiate\nas early as E9.0, prior to our sampling window (E10.5–E14.5), while late-stage myogenic cells at\nE15.5 form multinucleated myofibers that are challenging to isolate and analyze using sci-\nmtChIL-seq. These biological and technical factors likely contributed to the inconsistency\nbetween Lineage 1 and Lineage 2.\nIn light of these limitations, we instead focused the revised analysis on a MyoD-inducible\nNIH3T3 model, which provides a well-characterized and synchronized system for investigating\nMyoD- and Myog-dependent transcriptional dynamics. This model allowed us to resolve binding\ntransitions with greater clarity and reproducibility.\nWe believe this revision resolves the concern and improves the interpretability and robustness\nof the transcription factor binding analysis.\nOther comments:\nThere are numerous sections of the text that require major editing for clarity, accuracy, and to\nprovide enough information. Below are a few examples:\nFig 1i left RNAPIII\nFig 4f RNAII\nExtended Data Fig 3a left tH3K4me3\nWe thank the reviewer for identifying these labeling inaccuracies. In the revised manuscript, the\nfigures containing these annotations (original Fig. 1i, Fig. 4f, and Extended Data Fig. 3a) have\nbeen removed or substantially revised, and the specific panels referenced in the comment are no\nlonger included.\nWhile these labeling issues no longer persist in the current figures, we have reviewed all figure\nannotations throughout the manuscript to ensure consistency and accuracy. We greatly appreciate\nthe reviewer’s attention to detail, which contributed to improving the overall clarity of our\npresentation.\nReviewer #2 (Remarks to the Author):\nHarada and coauthors present multi-target ChIL-seq (mtChIL-seq) to simultaneously analyze the\nbinding of TFs, RNAPII, and other molecules to chromatins and the high-throughput extension\nof this chemistry to develop single-cell combinatorial indexing mtChIL-seq (sci-mtChIL-seq) that\ncan indirectly trace cell lineages with RNAPII and evaluate epigenomic status. They performed\nsci-mtChIL-seq on the lineage of embryonic myogenesis of mouse limbs to classify cell types\nbased on the binding state of RNAPII and simultaneously profile changes in the epigenomic state.\nThey also analyzed the role of TFs in cell lineage determination by invetigating MyoD-dependent\nepigenomic state development using in vitro models. Overall, this manuscripts presents some\nquite interesting data. The development of sci-mtChIL-seq itself in this work is largely based on\nthe ChIL–seq method developed by the same group, and provides a single-cell ChIL-seq profiling\nsolution with the split-pool barcoding strategy. There are some technical comments.\nThe authors provide solutions to both single-target sci-ChIL-seq and the multi-target sci-mtChIL-\nseq. But the first section of the manuscript only covers the single-target sci-ChIL-seq. Neither the\ncontext or the figures cover any details regarding how multiple targets can be profiled at the same\ntime. Did the authors incubate primary antibodies that target multiple antigens at the same time?\nIn the multi-target method (for example RNAPII and H3K27me3), are the sequencing reads\naggregates of RNAPII and histone signals, or the authors can stratify/demultiplex the RNAPII\nand H3K27me3 profiles?\nWe thank the reviewer for highlighting the need to clarify the structure and experimental design\nof our multi-target approach. In the revised manuscript, we have reorganized Fig. 1 to present\nboth the single-target and multi-target versions of the sci-ChIL-seq framework. To directly\naddress the reviewer’s questions:\n1. Does the first section cover multi-target profiling?\nYes, the revised Fig. 1e–i explicitly presents the design, validation, and output of multi-target sci-\nmtChIL-seq, including co-profiling of RNAPII and H3K4me3 in the same cells.\n Fig. 1b–d shows results from single-target sci-ChIL-seq.\n Fig. 1e–i introduce the multi-target method and demonstrates its capacity for target-specific\nread separation\n2. Were antibodies incubated together, and how was target identity preserved?\nIn the sci-mtChIL-seq protocol, we simultaneously incubate primary antibodies raised in\ndifferent species (e.g., mouse anti-RNAPII and rabbit anti-H3K4me3). These are then detected\nusing species-specific sci-ChIL probes, each loaded with distinct T5 barcode sequences. This\nsetup enables simultaneous detection without cross-reactivity and preserves target identity during\nbarcoding.\n3. Can the sequencing reads be stratified by target?\nYes. Because the T5 barcodes used in the probes are target-specific, the resulting sequencing reads\ncan be demultiplexed according to antibody identity. Therefore, RNAPII and histone modification\nsignals are not aggregated, but can be analyzed separately and independently even within the same\ncells.\nThese features are demonstrated in Fig. 1h–i, where we show the separation and target-specific\nread mapping across the genome. Additional clarifications have been added to the Results and\nMethods sections to describe this process in more detail.\nWe appreciate the reviewer’s helpful suggestion, which prompted us to improve the structure\nand clarity of this key methodological advance.\nAfter the bulk ChIL–seq was described in 2018, methods like CUT&RUN and CUT&Tag turn to\nbe the more accurate and robust for histone or epigenetic factor binding profiling. Can the authors\ncompare the data throughput, sensitivity, specificity, and performance with these related methods?\nWe appreciate the reviewer’s suggestion to compare sci-mtChIL-seq with related chromatin\nprofiling methods such as CUT&RUN and CUT&Tag. In the revised manuscript, we addressed\nthis point by focusing on data throughput, a broadly comparable and biologically context-\nindependent metric. Specifically, we show that sci-mtChIL-seq achieves comparable or higher\nnumbers of unique fragments per cell than recent multi-target methods including NTT-seq,\nMulTI-Tag, and uCoTarget (Revised manuscript, Page 5, lines 1–4).\nAdditionally, we performed internal validations of signal quality and reproducibility within sci-\nmtChIL-seq by comparing single- and multi-target formats (Fig. 1f–i), confirming consistent and\ntarget-specific signal profiles.\nWhile we agree that benchmarking against CUT&RUN and CUT&Tag in terms of sensitivity\nand specificity would be valuable, we did not perform direct comparisons for these metrics in this\nstudy. This is because such evaluations are highly dependent on the cell type, chromatin context,\nand antibody, and our study involved NIH3T3 and embryonic tissues, which differ from the\nstandardized systems (e.g., K562 or HEK293) commonly used in CUT&RUN/CUT&Tag\nbenchmarks. Without matched experimental setups, comparisons across methods may not yield\nmeaningful or interpretable results.\nMore details of this protocol should be described, such as why the authors chose to distribute a\ntotal of 25 cells per well in split-pool barcoding. Is FACS sorting required for this procedure? It\nis recommended to mention key procedures such as cell fixation, permeabilization and blocking\nin the main text to help readers get a better understanding of this technology.\nWe thank the reviewer for this valuable suggestion to improve the clarity and reproducibility of\nour protocol. In the revised manuscript, we have modified the main text to explicitly mention the\nkey preparatory steps of the sci-ChIL-seq workflow. Specifically, we added the following\nsentence to the description of Fig. 1a:\n“Prior to this, cells are fixed, permeabilized, and blocked following our previously established\nChIL-seq protocols[13,14].”\nThese steps correspond to 1% PFA fixation for 5 minutes, 1% Triton X-100 permeabilization for\n20 minutes, and blocking with Blocking One-P for 20 minutes, as detailed in the Methods section\n(Methods, Page 14 lines 11–21).\nIn addition, we clarified the rationale for using 25 cells per well during split-pool barcoding by\nadding:\n“We distributed 25 cells per well during split-pool barcoding to reduce the probability that two\ncells with the same T5 index are sorted into the same well. This approach minimizes index\ncollisions and doublet events, and is a standard practice in sci-based methods to preserve single-\ncell resolution.”\nLastly, we confirm that FACS sorting was used to isolate Hoechst-positive nuclei prior to\nbarcoding (Methods, Page 15 lines 10–14).\nWe appreciate the reviewer’s input, which helped us enhance the transparency and accessibility\nof the protocol in both the main text and the Methods section.\n[13] Harada, A. et al. A chromatin integration labelling method enables epigenomic profiling with\nlower input. Nat Cell Biol 21, (2019).\n[14] Handa, T. et al. Chromatin integration labeling for mapping DNA-binding proteins and\nmodifications with low input. Nat Protoc 15, (2020).\nThe authors mention they use homemade Tn5. Please describe the Tn5 preparation in\nmethod and add relevant citations for the protocol.\nWe thank the reviewer for this pointing this out. In the revised manuscript, we have added a\ndescription of the preparation and use of the in-house Tn5 transposase in the Methods section.\nSpecifically, we now state that the enzyme was prepared in-house and previously published\nprotocol [1], and was loaded with mosaic-end sequences prior to tagmentation.\nThis addition improves the transparency and reproducibility of our experimental workflow. The\nrevised text reads:\n“The Tn5 transposase was produced in-house following a previously described protocol [41] and\nwas pre-loaded with mosaic end sequences.”\nWe also note that the production of the Tn5 enzyme was performed by the authors (S.S. and\nH.Ku.), as stated in the Author Contributions section.\n[41] Sato, S. et al. Biochemical analysis of nucleosome targeting by Tn5 transposase. Open Biol.\n9, 190116 (2019).\nReviewer #3 (Remarks to the Author):\nHarada et al develop a new assay to probe RNA Pol II binding and epigenomic factors in single\ncells. They then use this tool to evaluate binding dynamics during myogenesis within the\ndeveloping murine limb. While this is an exciting tool, it’s not clear what new information is\nincluded in this manuscript. The results are not clearly articulated and not discussed in the context\nof what is known.\nExamples of text that are unclear or that require context:\n-The text accommodating Fig. 1i-j says, “..association between cell populations represented by\nRNAPII binding and histone modifications during limb development.” But what exactly does that\ntell you? Also, “..the cell population represented by H3K27me3 was..” How do you define large?\nAnd why is that significant? These are just two examples of a common theme from the manuscript.\nWe thank the reviewer for these insightful questions. In the revised manuscript, we have\naddressed these concerns by removing both original Fig. 1i and 1j, and reorganizing the associated\nanalysis to improve biological interpretability and eliminate ambiguous language. The relevant\nco-profiling content is now presented in Supplementary Fig. 3, using RNAPII-based integration\nacross datasets.\nRegarding the first part of the comment, the phrase “association between cell populations” in\nthe original Fig. 1i–j referred to an analysis of whether chromatin states defined by histone\nmodifications (H3K4me3 or H3K27me3) correspond to transcriptional states inferred from\nRNAPII profiles in the same cells. To visualize this, we used UMAP for each epigenomic factor\nseparately and connected cells across modalities, enabling us to assess inter-modality consistency.\nIn the RNAPII–H3K4me3 dataset (original Fig. 1i), we observed strong correspondence\nbetween modalities, indicating that transcriptionally active cells also carry H3K4me3 signals—\nconsistent with known biology. In contrast, in the RNAPII–H3K27me3 dataset (original Fig. 1j),\nthe H3K27me3 signal was lacked RNAPII-defined cell state specificity. We originally described\nthis as “a large cell population,” but we now recognize that this term was ambiguous and possibly\nmisleading. As a result, we have removed Fig. 1j and the associated wording.\nTo resolve this and improve interpretability, we conducted multi-target integration using RNAPII\nas a common axis, allowing comparison of H3K4me3 and H3K27me3 profiles based on common\nspace. As shown in Supplementary Fig. 3, this approach reveals that while active and repressive\nchromatin marks in vary cell state, the transcriptional backbone (RNAPII) provides a stable and\ninterpretable structure for defining cell states across datasets.\nWe also recognize the reviewer’s broader point that these issues reflect a general theme of\nimprecise terminology or unclear interpretation throughout the manuscript. In response, we have\ncarefully reviewed and revised other parts of the text to improve consistency, clarify the biological\nmeaning of our terms, and avoid ambiguous or overstated conclusions.\nWe appreciate the reviewer’s comment, which helped us refine both the conceptual clarity and\npresentation of our results.\n-Regarding the MyoD binding analyses, what does Fig. 3D show? The text first says a motif\nenrichment was performed but then it says that the ‘myocyte lineage had the greater deviation\nscore.” What is a deviation score?\nWe thank the reviewer for this thoughtful comment. In the revised manuscript, we have removed\nthe original Fig. 3d and 3g, which presented trajectory-based lineage inference and transcription\nfactor motif activity estimation, respectively. These analyses were excluded to improve clarity\nand to focus the manuscript on experimentally supported results.\nNevertheless, we would like to clarify the original intent and interpretation of these analyses.\nThe goal was to infer the functional activity of MyoD during myogenesis by evaluating the\nenrichment of MyoD-binding motifs across different pseudotime-defined cell clusters.\nTo quantify the enrichment levels of TF motifs around the MyoD bound sites, we applied\nchromVAR, which is developed for ATAC-seq data to calculate a deviation score reflecting how\nmuch motif accessibility deviates from a background model across datasets. In this way, the\ndeviation score serves as a quantitative proxy for transcription factor activity, allowing relative\ncomparisons of motif involvement across datasets.\nIn our case, MyoD bound sites were used as chromVAR inpu. A higher deviation score in the\nmyocytes indicated that these regions showed greater variability and were more specifically\nengaged in those cells—consistent with MyoD's known role in muscle differentiation. The use of\nchromVAR and deviation scores in this context is well established in single-cell chromatin\nprofiling, including scATAC-seq and ChIL-based methods.\nThat said, we agree with the reviewer that motif-level inference is not equivalent to direct\nbinding evidence, and that interpretation based on motif analysis should be supported by\nfunctional validation. To avoid potential overinterpretation and streamline the manuscript, we\ndecided to omit this analysis and its associated figures from the revised version.\n-The authors also reference the term lineage tracing multiple times but as described, their assay\ndoes not evaluate lineage as it is traditionally defined. Trajectory inference allows you to infer\nfate but it does not establish lineage identity.\nWe appreciate the reviewer’s attention to the precision of terminology. In the revised manuscript,\nwe carefully avoided the use of the term lineage tracing, and instead described our approach using\ntrajectory inference or epigenomic trajectory to infer temporal regulatory order or differentiation\ndirectionality based on chromatin states. Our method does not aim to trace heritable cell lineages\nbut rather captures state transitions inferred from chromatin dynamics. We agree that “lineage\ntracing” is reserved for approaches that follow cell progeny using genetic labeling, and we have\nensured our revised text reflects this conceptual distinction."}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors have made a substantial revision in response to my previous comments. The revised manuscript now clearly explains the principles and implementation of sci-mtChIL-seq, the strategy for multiplexed antibody labeling and barcode-based read separation, which was not well detailed in the original version. The updated figures and restructured text improve clarity, and key methodological steps such as fixation, barcoding, and rationale for cell loading are now described in details. This reviewer also appreciates the addition of a public analysis pipeline and the efforts to improve figure labeling and consistency. While a direct performance comparison to related methods such as CUT&RUN and CUT&Tag remains limited, the authors provide reasonable justification and include throughput benchmarking against other multi-target platforms. Overall, the revised manuscript much better articulates the value of sci-mtChIL-seq as a flexible and scalable tool for profiling chromatin states and transcriptional regulation at single-cell level.", "reviewer_4": "The Reviewer #1's comments are addressed in this revision.", "reviewer_5": "Reviewers' comments:\nReviewer #1 (Remarks to the Author):\nThis study by Harada et al. investigates the regulatory mechanisms underlying myogenic\ndifferentiation during mouse embryonic development using single-cell multi-omics techniques.\nThe authors profiled chromatin interactions and transcriptional activities in myogenic progenitor\ncells across key developmental stages. They employed a sci-mtChIL-seq methodologies to map\nchromatin accessibility, histone modifications, and transcription factor binding dynamics with\nhigh spatial and temporal resolution. Their findings reveal dynamic changes in chromatin\narchitecture and transcriptional regulation orchestrated by key regulatory factors such as MyoD\nand Myog. Moreover, the study highlights the utility of single-cell multi-omics approaches in\nunraveling complex gene regulatory networks during embryonic development.\nSpecific comments are as follows:\n1. All conclusions drawn in the manuscript were not validated experimentally, which significantly\nundermines the strength of their claims.\nWe thank the reviewer for this important point. In the revised manuscript, we addressed this\nconcern by demonstrating that our analytical framework is capable of recapitulating well-\nestablished, experimentally validated regulatory phenomena. Specifically:\n In Fig. 3, our data successfully reproduce the known sequential binding of MyoD followed\nby Myog during myogenic differentiation, consistent with prior perturbation-based studies\n(Cao et al., EMBO J. 2006, PMID: 16437161; Ohkawa et al., EMBO J. 2006, PMID:\n16424906).\n In Fig. 4, our method reconstructs the temporal order of chromatin regulator recruitment at\nmuscle gene loci, in line with established findings on regulatory coordination in muscle genes\n(Hernández-Hernández et al., Front. Genet. 2020, PMID: 33193700; Blum et al., Genes Dev.\n2012, PMID: 23249738).\nThese validations demonstrate that our framework can recover biologically meaningful\nregulatory cascades at single-cell resolution, even in the absence of new perturbation experiments.\nIn the original manuscript, we had focused on demonstrating the utility of our method through\ntechnically consistent patterns derived from known regulatory relationships. However, we now\nrecognize that certain interpretations lacked sufficient biological context or direct support. As a\nresult, we have substantially revised the manuscript, removing under-supported analyses and\nincorporating additional layers of evidence where appropriate.\nFor example, in Supplementary Fig. 5e–f, we show that at the Mef2a locus, the earliest\naccumulation of chromatin modifiers (Brg1 and H3K27ac) occurs at a distal enhancer annotated\nin the ENCODE database, prior to MyoD binding at the promoter, H3K4me3 deposition, Myog\nbinding, and RNAPII elongation. This sequential pattern suggests a temporally ordered cascade\nof regulatory events.\nFurthermore, bulk ATAC-seq data confirm that this distal enhancer region is already accessible\nat the undifferentiated (MyoD–) stage, suggesting transcriptional priming through pre-established\nchromatin accessibility.\nTogether, these findings support a model in which MyoD sequentially engages enhancer and\npromoter regions to promote active chromatin formation, followed by Myog and RNAPII\nrecruitment. The reproducibility of these patterns, their consistency with known enhancer biology,\nand alignment with orthogonal data sources provide strong technical and biological validation of\nour conclusions.\n2. The method description was not informative. This is critical because the conclusion could lead\nto overfitting.\nWe agree with the reviewer that insufficient methodological clarity may weaken the\ninterpretation of our results, particularly with regard to the risk of overfitting. In the revised\nmanuscript, we have clarified the rationale and implementation of trajectory-based visualizations,\nremoved figures that were difficult to interpret due to dataset-specific pseudotime assignments,\nand made the complete analysis pipeline publicly available to ensure transparency and\nreproducibility.\nSpecifically, we inferd that the reviewer’s concern may relate to the trend estimation in the\noriginal Fig. 3e,f and Fig. 4c,g, where we plotted changes in factor accumulation over pseudotime\nusing regression curves. In particular, Fig. 4c,g showed accumulation levels averaged within\npseudotime windows to reduce noise and enhance robustness.\nTo mitigate overfitting, we used LOESS (locally estimated scatterplot smoothing) with a large\nspan parameter, which limits model flexibility. While this reduces temporal resolution, it helps\nprevent the regression from capturing spurious fluctuations. Thus, we believe the risk of\noverfitting in these visualizations was minimal.\nNevertheless, we acknowledge that interpreting trends was complicated by the use of\nindependently computed pseudotime for each dataset, which limited direct comparisons. To\naddress this, we have removed the original panels and, in the revised Fig. 4, present an alternative\nanalysis that integrates multiple datasets using RNAPII as a shared reference. This analysis\nevaluates factor accumulation at individual gene loci, offering a more intuitive and interpretable\nview of regulatory order without relying on pseudotime inference or regression-based smoothing.\nFurthermore, we have made the full analysis pipeline—including data preprocessing,\npseudotime inference, GAM smoothing parameters, and visualization scripts—publicly available\nvia GitHub (see Code availability section). While the revised text outlines the core methodology,\nadditional technical details are available in the repository to ensure full transparency.\nWe thank the reviewer for highlighting this important issue and prompting us to refine the\npresentation and interpretation of our trajectory-based analyses.\n3. Page 4 line 4 It is not clear how to set cut-off values and cell-derived indices.\nWe agree with the reviewer that the method for defining cut-off thresholds and identifying cell-\nderived indices should be clearly explained. In the revised manuscript, we now describe how these\nthresholds were estimated in a data-driven and dataset-specific manner using a statistical model.\nSpecifically, for barcode demultiplexing, we used a combinatorial split-pool indexing strategy,\nin which unique cell identities were assigned via multiple rounds of barcoding. To remove\nbarcodes with low read counts likely derived from ambient background or non-cell-containing\nindices, we used a two-component Gaussian mixture model and estimated the threshold values as\nthe intersection point between the two distributions. Additionally, we inferred doublet probability\nbased on the frequency of occurrence per index and exclude potential doublets or multiplets.\n4. Page 4 line 8 “The genomic localization patterns of RNAPII from RNAPII-H3K4me3 and\nRNAPII-H3K27me3 were roughly similar between stages.” How can this result be explained, and\ndoes it indicate that the localization of PolII is not related to the levels of H3K4me3 and\nH3K27me3?\nWe thank the reviewer for pointing out the ambiguity in our original statement. In response, we\nhave removed the sentence in question from the revised manuscript to avoid misinterpretation.\nIn the revised version, we instead emphasize that RNAPII signals from both RNAPII-H3K4me3\nand RNAPII-H3K27me3 datasets were successfully integrated into a common space using\nHarmony, and showed consistent clustering without major batch-specific segregation\n(Supplementary Fig. 3b). This supports the conclusion that RNAPII localization reflects shared\ntranscriptional activity across stages, while also being compatible with the chromatin context (e.g.,\nactive or poised states).\nWe hope this clarification resolves the concern and makes our interpretation more precise.\n5. Page 5 line 4 What is stage-dependent bias?\nWe thank the reviewer for this important question. In the original manuscript, the term “stage-\ndependent bias” was used to describe an observation in the H3K27me3 dataset, where cell\nclustering patterns corresponded more strongly to developmental stage than to cell type identity.\nHowever, we now recognize that this terminology may have been ambiguous and potentially\nmisleading.\nSpecifically, in the H3K27me3 dataset derived from embryonic limb tissues, UMAP\nvisualization showed the cluster structure that aligned closely with embryonic day (e.g., E10.5 vs.\nE12.5), rather than with transcriptionally distinct cell populations. This pattern likely arises from\nbroad H3K27me3-marked heterochromatin domains, which can be partially shared across\nlineages but shift globally during developmental progression. Thus, the stage-dependent pattern\nreflects a biological phenomenon related to chromatin state transitions, not a technical artifact.\nHowever, we agree with the reviewer that the term “bias” may unintentionally suggest a\nconfounding or technical error. To avoid this confusion, we have removed the term and its\nassociated interpretation from the revised manuscript. This decision was made not due to any data\nquality issue, but to improve clarity and avoid misinterpretation. We emphasize that the embryo\ndatasets remain a core component of our analysis, particularly for RNAPII-based cell state\nclassification and trajectory inference. The removal of this terminology does not reflect a\ndismissal of the embryo data, but rather a refinement of how we present their biological\ninterpretation.\nWe thank the reviewer again for drawing our attention to the potential ambiguity in our wording.\n6. Fig 3e-f MyoD are inconsistent between Lineage1 and Lineage2.\nWe thank the reviewer for this valuable comment. We have addressed this issue by removing the\noriginal MyoD binding analysis in embryonic samples from the revised manuscript.\nAs we noted in the first submission, the number of MyoD-expressing cells captured from\nembryonic tissues was limited, making reliable trajectory reconstruction and quantitative\ncomparison across lineages technically difficult. Moreover, MyoD expression is known to initiate\nas early as E9.0, prior to our sampling window (E10.5–E14.5), while late-stage myogenic cells at\nE15.5 form multinucleated myofibers that are challenging to isolate and analyze using sci-\nmtChIL-seq. These biological and technical factors likely contributed to the inconsistency\nbetween Lineage 1 and Lineage 2.\nIn light of these limitations, we instead focused the revised analysis on a MyoD-inducible\nNIH3T3 model, which provides a well-characterized and synchronized system for investigating\nMyoD- and Myog-dependent transcriptional dynamics. This model allowed us to resolve binding\ntransitions with greater clarity and reproducibility.\nWe believe this revision resolves the concern and improves the interpretability and robustness\nof the transcription factor binding analysis.\nOther comments:\nThere are numerous sections of the text that require major editing for clarity, accuracy, and to\nprovide enough information. Below are a few examples:\nFig 1i left RNAPIII\nFig 4f RNAII\nExtended Data Fig 3a left tH3K4me3\nWe thank the reviewer for identifying these labeling inaccuracies. In the revised manuscript, the\nfigures containing these annotations (original Fig. 1i, Fig. 4f, and Extended Data Fig. 3a) have\nbeen removed or substantially revised, and the specific panels referenced in the comment are no\nlonger included.\nWhile these labeling issues no longer persist in the current figures, we have reviewed all figure\nannotations throughout the manuscript to ensure consistency and accuracy. We greatly appreciate\nthe reviewer’s attention to detail, which contributed to improving the overall clarity of our\npresentation.\nReviewer #2 (Remarks to the Author):\nHarada and coauthors present multi-target ChIL-seq (mtChIL-seq) to simultaneously analyze the\nbinding of TFs, RNAPII, and other molecules to chromatins and the high-throughput extension\nof this chemistry to develop single-cell combinatorial indexing mtChIL-seq (sci-mtChIL-seq) that\ncan indirectly trace cell lineages with RNAPII and evaluate epigenomic status. They performed\nsci-mtChIL-seq on the lineage of embryonic myogenesis of mouse limbs to classify cell types\nbased on the binding state of RNAPII and simultaneously profile changes in the epigenomic state.\nThey also analyzed the role of TFs in cell lineage determination by invetigating MyoD-dependent\nepigenomic state development using in vitro models. Overall, this manuscripts presents some\nquite interesting data. The development of sci-mtChIL-seq itself in this work is largely based on\nthe ChIL–seq method developed by the same group, and provides a single-cell ChIL-seq profiling\nsolution with the split-pool barcoding strategy. There are some technical comments.\nThe authors provide solutions to both single-target sci-ChIL-seq and the multi-target sci-mtChIL-\nseq. But the first section of the manuscript only covers the single-target sci-ChIL-seq. Neither the\ncontext or the figures cover any details regarding how multiple targets can be profiled at the same\ntime. Did the authors incubate primary antibodies that target multiple antigens at the same time?\nIn the multi-target method (for example RNAPII and H3K27me3), are the sequencing reads\naggregates of RNAPII and histone signals, or the authors can stratify/demultiplex the RNAPII\nand H3K27me3 profiles?\nWe thank the reviewer for highlighting the need to clarify the structure and experimental design\nof our multi-target approach. In the revised manuscript, we have reorganized Fig. 1 to present\nboth the single-target and multi-target versions of the sci-ChIL-seq framework. To directly\naddress the reviewer’s questions:\n1. Does the first section cover multi-target profiling?\nYes, the revised Fig. 1e–i explicitly presents the design, validation, and output of multi-target sci-\nmtChIL-seq, including co-profiling of RNAPII and H3K4me3 in the same cells.\n Fig. 1b–d shows results from single-target sci-ChIL-seq.\n Fig. 1e–i introduce the multi-target method and demonstrates its capacity for target-specific\nread separation\n2. Were antibodies incubated together, and how was target identity preserved?\nIn the sci-mtChIL-seq protocol, we simultaneously incubate primary antibodies raised in\ndifferent species (e.g., mouse anti-RNAPII and rabbit anti-H3K4me3). These are then detected\nusing species-specific sci-ChIL probes, each loaded with distinct T5 barcode sequences. This\nsetup enables simultaneous detection without cross-reactivity and preserves target identity during\nbarcoding.\n3. Can the sequencing reads be stratified by target?\nYes. Because the T5 barcodes used in the probes are target-specific, the resulting sequencing reads\ncan be demultiplexed according to antibody identity. Therefore, RNAPII and histone modification\nsignals are not aggregated, but can be analyzed separately and independently even within the same\ncells.\nThese features are demonstrated in Fig. 1h–i, where we show the separation and target-specific\nread mapping across the genome. Additional clarifications have been added to the Results and\nMethods sections to describe this process in more detail.\nWe appreciate the reviewer’s helpful suggestion, which prompted us to improve the structure\nand clarity of this key methodological advance.\nAfter the bulk ChIL–seq was described in 2018, methods like CUT&RUN and CUT&Tag turn to\nbe the more accurate and robust for histone or epigenetic factor binding profiling. Can the authors\ncompare the data throughput, sensitivity, specificity, and performance with these related methods?\nWe appreciate the reviewer’s suggestion to compare sci-mtChIL-seq with related chromatin\nprofiling methods such as CUT&RUN and CUT&Tag. In the revised manuscript, we addressed\nthis point by focusing on data throughput, a broadly comparable and biologically context-\nindependent metric. Specifically, we show that sci-mtChIL-seq achieves comparable or higher\nnumbers of unique fragments per cell than recent multi-target methods including NTT-seq,\nMulTI-Tag, and uCoTarget (Revised manuscript, Page 5, lines 1–4).\nAdditionally, we performed internal validations of signal quality and reproducibility within sci-\nmtChIL-seq by comparing single- and multi-target formats (Fig. 1f–i), confirming consistent and\ntarget-specific signal profiles.\nWhile we agree that benchmarking against CUT&RUN and CUT&Tag in terms of sensitivity\nand specificity would be valuable, we did not perform direct comparisons for these metrics in this\nstudy. This is because such evaluations are highly dependent on the cell type, chromatin context,\nand antibody, and our study involved NIH3T3 and embryonic tissues, which differ from the\nstandardized systems (e.g., K562 or HEK293) commonly used in CUT&RUN/CUT&Tag\nbenchmarks. Without matched experimental setups, comparisons across methods may not yield\nmeaningful or interpretable results.\nMore details of this protocol should be described, such as why the authors chose to distribute a\ntotal of 25 cells per well in split-pool barcoding. Is FACS sorting required for this procedure? It\nis recommended to mention key procedures such as cell fixation, permeabilization and blocking\nin the main text to help readers get a better understanding of this technology.\nWe thank the reviewer for this valuable suggestion to improve the clarity and reproducibility of\nour protocol. In the revised manuscript, we have modified the main text to explicitly mention the\nkey preparatory steps of the sci-ChIL-seq workflow. Specifically, we added the following\nsentence to the description of Fig. 1a:\n“Prior to this, cells are fixed, permeabilized, and blocked following our previously established\nChIL-seq protocols[13,14].”\nThese steps correspond to 1% PFA fixation for 5 minutes, 1% Triton X-100 permeabilization for\n20 minutes, and blocking with Blocking One-P for 20 minutes, as detailed in the Methods section\n(Methods, Page 14 lines 11–21).\nIn addition, we clarified the rationale for using 25 cells per well during split-pool barcoding by\nadding:\n“We distributed 25 cells per well during split-pool barcoding to reduce the probability that two\ncells with the same T5 index are sorted into the same well. This approach minimizes index\ncollisions and doublet events, and is a standard practice in sci-based methods to preserve single-\ncell resolution.”\nLastly, we confirm that FACS sorting was used to isolate Hoechst-positive nuclei prior to\nbarcoding (Methods, Page 15 lines 10–14).\nWe appreciate the reviewer’s input, which helped us enhance the transparency and accessibility\nof the protocol in both the main text and the Methods section.\n[13] Harada, A. et al. A chromatin integration labelling method enables epigenomic profiling with\nlower input. Nat Cell Biol 21, (2019).\n[14] Handa, T. et al. Chromatin integration labeling for mapping DNA-binding proteins and\nmodifications with low input. Nat Protoc 15, (2020).\nThe authors mention they use homemade Tn5. Please describe the Tn5 preparation in\nmethod and add relevant citations for the protocol.\nWe thank the reviewer for pointing this out. In the revised manuscript, we have added a\ndescription of the preparation and use of the in-house Tn5 transposase in the Methods section.\nSpecifically, we now state that the enzyme was prepared in-house and previously published\nprotocol [1], and was loaded with mosaic-end sequences prior to tagmentation.\nThis addition improves the transparency and reproducibility of our experimental workflow. The\nrevised text reads:\n“The Tn5 transposase was produced in-house following a previously described protocol [41] and\nwas pre-loaded with mosaic end sequences.”\nWe also note that the production of the Tn5 enzyme was performed by the authors (S.S. and\nH.Ku.), as stated in the Author Contributions section.\n[41] Sato, S. et al. Biochemical analysis of nucleosome targeting by Tn5 transposase. Open Biol.\n9, 190116 (2019).\nReviewer #3 (Remarks to the Author):\nHarada et al develop a new assay to probe RNA Pol II binding and epigenomic factors in single\ncells. They then use this tool to evaluate binding dynamics during myogenesis within the\ndeveloping murine limb. While this is an exciting tool, it’s not clear what new information is\nincluded in this manuscript. The results are not clearly articulated and not discussed in the context\nof what is known.\nExamples of text that are unclear or that require context:\n-The text accommodating Fig. 1i-j says, “..association between cell populations represented by\nRNAPII binding and histone modifications during limb development.” But what exactly does that\ntell you? Also, “..the cell population represented by H3K27me3 was..” How do you define large?\nAnd why is that significant? These are just two examples of a common theme from the manuscript.\nWe thank the reviewer for these insightful questions. In the revised manuscript, we have\naddressed these concerns by removing both original Fig. 1i and 1j, and reorganizing the associated\nanalysis to improve biological interpretability and eliminate ambiguous language. The relevant\nco-profiling content is now presented in Supplementary Fig. 3, using RNAPII-based integration\nacross datasets.\nRegarding the first part of the comment, the phrase “association between cell populations” in\nthe original Fig. 1i–j referred to an analysis of whether chromatin states defined by histone\nmodifications (H3K4me3 or H3K27me3) correspond to transcriptional states inferred from\nRNAPII profiles in the same cells. To visualize this, we used UMAP for each epigenomic factor\nseparately and connected cells across modalities, enabling us to assess inter-modality consistency.\nIn the RNAPII–H3K4me3 dataset (original Fig. 1i), we observed strong correspondence\nbetween modalities, indicating that transcriptionally active cells also carry H3K4me3 signals—\nconsistent with known biology. In contrast, in the RNAPII–H3K27me3 dataset (original Fig. 1j),\nthe H3K27me3 signal was lacked RNAPII-defined cell state specificity. We originally described\nthis as “a large cell population,” but we now recognize that this term was ambiguous and possibly\nmisleading. As a result, we have removed Fig. 1j and the associated wording.\nTo resolve this and improve interpretability, we conducted multi-target integration using RNAPII\nas a common axis, allowing comparison of H3K4me3 and H3K27me3 profiles based on common\nspace. As shown in Supplementary Fig. 3, this approach reveals that while active and repressive\nchromatin marks in vary cell state, the transcriptional backbone (RNAPII) provides a stable and\ninterpretable structure for defining cell states across datasets.\nWe also recognize the reviewer’s broader point that these issues reflect a general theme of\nimprecise terminology or unclear interpretation throughout the manuscript. In response, we have\ncarefully reviewed and revised other parts of the text to improve consistency, clarify the biological\nmeaning of our terms, and avoid ambiguous or overstated conclusions.\nWe appreciate the reviewer’s comment, which helped us refine both the conceptual clarity and\npresentation of our results.\n-Regarding the MyoD binding analyses, what does Fig. 3D show? The text first says a motif\nenrichment was performed but then it says that the ‘myocyte lineage had the greater deviation\nscore.” What is a deviation score?\nWe thank the reviewer for this thoughtful comment. In the revised manuscript, we have\nremoved the original Fig. 3d and 3g, which presented trajectory-based lineage inference and transcription\nfactor motif activity estimation, respectively. These analyses were excluded to improve clarity\nand to focus the manuscript on experimentally supported results.\nNevertheless, we would like to clarify the original intent and interpretation of these analyses.\nThe goal was to infer the functional activity of MyoD during myogenesis by evaluating the\nenrichment of MyoD-binding motifs across different pseudotime-defined cell clusters.\nTo quantify the enrichment levels of TF motifs around the MyoD bound sites, we applied\nchromVAR, which is developed for ATAC-seq data to calculate a deviation score reflecting how\nmuch motif accessibility deviates from a background model across datasets. In this way, the\ndeviation score serves as a quantitative proxy for transcription factor activity, allowing relative\ncomparisons of motif involvement across datasets.\nIn our case, MyoD bound sites were used as chromVAR inpu. A higher deviation score in the\nmyocytes indicated that these regions showed greater variability and were more specifically\nengaged in those cells—consistent with MyoD's known role in muscle differentiation. The use of\nchromVAR and deviation scores in this context is well established in single-cell chromatin\nprofiling, including scATAC-seq and ChIL-based methods.\nThat said, we agree with the reviewer that motif-level inference is not equivalent to direct\nbinding evidence, and that interpretation based on motif analysis should be supported by\nfunctional validation. To avoid potential overinterpretation and streamline the manuscript, we\ndecided to omit this analysis and its associated figures from the revised version.\n-The authors also reference the term lineage tracing multiple times but as described, their assay\ndoes not evaluate lineage as it is traditionally defined. Trajectory inference allows you to infer\nfate but it does not establish lineage identity.\nWe appreciate the reviewer’s attention to the precision of terminology. In the revised manuscript,\nwe carefully avoided the use of the term lineage tracing, and instead described our approach using\ntrajectory inference or epigenomic trajectory to infer temporal regulatory order or differentiation\ndirectionality based on chromatin states. Our method does not aim to trace heritable cell lineages\nbut rather captures state transitions inferred from chromatin dynamics. We agree that “lineage\ntracing” is reserved for approaches that follow cell progeny using genetic labeling, and we have\nensured our revised text reflects this conceptual distinction."}, "author_response": ""}]

Decision letter