ULDNA: integrating unsupervised multi-source language models with LSTM-attention network for high-accuracy protein-DNA binding site prediction.

Efficient and accurate recognition of protein-DNA interactions is vital for understanding the molecular mechanisms of related biological processes and further guiding drug discovery. Although the current experimental protocols are the most precise way to determine protein-DNA binding sites, they tend to be labor-intensive and time-consuming. There is an immediate need to design efficient computational approaches for predicting DNA-binding sites. Here, we proposed ULDNA, a new deep-learning model, to deduce DNA-binding sites from protein sequences. This model leverages an LSTM-attention architecture, embedded with three unsupervised language models that are pre-trained on large-scale sequences from multiple database sources. To prove its effectiveness, ULDNA was tested on 229 protein chains with experimental annotation of DNA-binding sites. Results from computational experiments revealed that ULDNA significantly improves the accuracy of DNA-binding site prediction in comparison with 17 state-of-the-art methods. In-depth data analyses showed that the major strength of ULDNA stems from employing three transformer language models. Specifically, these language models capture complementary feature embeddings with evolution diversity, in which the complex DNA-binding patterns are buried. Meanwhile, the specially crafted LSTM-attention network effectively decodes evolution diversity-based embeddings as DNA-binding results at the residue level. Our findings demonstrated a new pipeline for predicting DNA-binding sites on a large scale with high accuracy from protein sequence alone.

Cooperative action of separate interaction domains promotes high-affinity DNA binding of Arabidopsis thaliana ARF transcription factors.

The signaling molecule auxin coordinates many growth and development processes in plants, mainly through modulating gene expression. Transcriptional response is mediated by the family of auxin response factors (ARF). Monomers of this family recognize a DNA motif and can homodimerize through their DNA-binding domain (DBD), enabling cooperative binding to an inverted binding site. Most ARFs further contain a C-terminal PB1 domain that is capable of homotypic interactions and mediating interactions with Aux/IAA repressors. Given the dual role of the PB1 domain, and the ability of both DBD and PB1 domain to mediate dimerization, a key question is how these domains contribute to DNA-binding specificity and affinity. So far, ARF-ARF and ARF-DNA interactions have mostly been approached using qualitative methods that do not provide a quantitative and dynamic view on the binding equilibria. Here, we utilize a DNA binding assay based on single-molecule Förster resonance energy transfer (smFRET) to study the affinity and kinetics of the interaction of several Arabidopsis thaliana ARFs with an IR7 auxin-responsive element (AuxRE). We show that both DBD and PB1 domains of AtARF2 contribute toward DNA binding, and we identify ARF dimer stability as a key parameter in defining binding affinity and kinetics across AtARFs. Lastly, we derived an analytical solution for a four-state cyclic model that explains both the kinetics and the affinity of the interaction between AtARF2 and IR7. Our work demonstrates that the affinity of ARFs toward composite DNA response elements is defined by dimerization equilibrium, identifying this as a key element in ARF-mediated transcriptional activity.

5-Formylcytosine mediated DNA-peptide cross-link induces predominantly semi-targeted mutations in both Escherichia coli and human cells.

Histone proteins can become trapped on DNA in the presence of 5-formylcytosine (5fC) to form toxic DNA-protein conjugates. Their repair may involve proteolytic digestion resulting in DNA-peptide cross-links (DpCs). Here, we have investigated replication of a model DpC comprised of an 11-mer peptide (NH2-GGGKGLGK∗GGA) containing an oxy-lysine residue (K∗) conjugated to 5fC in DNA. Both CXG and CXT (where X = 5fC-DpC) sequence contexts were examined. Replication of both constructs gave low viability (<10%) in Escherichia coli, whereas TLS efficiency was high (72%) in HEK 293T cells. In E. coli, the DpC was bypassed largely error-free, inducing only 2 to 3% mutations, which increased to 4 to 5% with SOS. For both sequences, semi-targeted mutations were dominant, and for CXG, the predominant mutations were G→T and G→C at the 3'-base to the 5fC-DpC. In HEK 293T cells, 7 to 9% mutations occurred, and the dominant mutations were the semi-targeted G → T for CXG and T → G for CXT. These mutations were reduced drastically in cells deficient in hPol η, hPol ι or hPol ζ, suggesting a role of these TLS polymerases in mutagenic TLS. Steady-state kinetics studies using hPol η confirmed that this polymerase induces G → T and T → G transversions at the base immediately 3' to the DpC. This study reveals a unique replication pattern of 5fC-conjugated DpCs, which are bypassed largely error-free in both E. coli and human cells and induce mostly semi-targeted mutations at the 3' position to the lesion.

Dynamic stem-loop extension by Pol θ and templated insertion during DNA repair.

Theta-mediated end joining (TMEJ) is critical for survival of cancer cells when other DNA double-stranded break repair pathways are impaired. Human DNA polymerase theta (Pol θ) can extend ssDNA oligonucleotides, but little is known about preferred substrates and mechanism. We show that Pol θ can extend both ssDNA and RNA substrates by unimolecular stem-loop synthesis initiated by only two 3' terminal base pairs. Given sufficient time, Pol θ uses alternative pairing configurations that greatly expand the repertoire of sequence outcomes. Further primer-template adjustments yield low-fidelity outcomes when the nucleotide pool is imbalanced. Unimolecular stem-loop synthesis competes with bimolecular end joining, even when a longer terminal microhomology for end joining is available. Both reactions are partially suppressed by the ssDNA-binding protein replication protein A. Protein-primer grasp residues that are specific to Pol θ are needed for rapid stem-loop synthesis. The ability to perform stem-loop synthesis from a minimally paired primer is rare among human DNA polymerases, but we show that human DNA polymerases Pol η and Pol λ can catalyze related reactions. Using purified human Pol θ, we reconstituted in vitro TMEJ incorporating an insertion arising from a stem-loop extension. These activities may help explain TMEJ repair events that include inverted repeat sequences.

A real-time biochemical assay for quantitative analyses of APOBEC-catalyzed DNA deamination.

Over the past decade, the connection between APOBEC3 cytosine deaminases and cancer mutagenesis has become increasingly apparent. This growing awareness has created a need for biochemical tools that can be used to identify and characterize potential inhibitors of this enzyme family. In response to this challenge, we have developed a Real-time APOBEC3-mediated DNA Deamination assay. This assay offers a single-step set-up and real-time fluorescent read-out, and it is capable of providing insights into enzyme kinetics. The assay also offers a high-sensitivity and easily scalable method for identifying APOBEC3 inhibitors. This assay serves as a crucial addition to the existing APOBEC3 biochemical and cellular toolkit and possesses the versatility to be readily adapted into a high-throughput format for inhibitor discovery.

The multifarious MerR family of transcriptional regulators.

The MerR family of transcriptional regulators includes a variety of bacterial cytoplasmic proteins that respond to a wide range of signals, including toxins, metal ions, and endogenous metabolites. Its best-characterized members share similar structural and functional features with the family founder, the mercury sensor MerR, although most of them do not respond to metal ions. The group of "canonical" MerR homologs displays common molecular mechanisms for controlling the transcriptional activation of their target genes in response to inducer signals. This includes the recognition of distinctive operator sequences located at suboptimal σ70 -dependent promoters. Interestingly, an increasing number of proteins assigned to the MerR family based on their DNA-binding domain do not match in structure, sequence, or mode of action with any of the canonical MerR-like regulators. Here, we analyzed several members of the family, including this last group. Based on a phylogenetic analysis, and similarities in structural/functional features and position of their target operators relative to the promoter elements, we propose to assign these "atypical/divergent" MerR regulators to a phylogenetically separated group. These atypical/divergent homologs represent a new class of transcriptional regulators with novel regulatory mechanisms.

Redefining the bacteriophage mv4 site-specific recombination system and the sequence specificity of its attB and core-attP sites.

Through their involvement in the integration and excision of a large number of mobile genetic elements, such as phages and integrative and conjugative elements (ICEs), site-specific recombination systems based on heterobivalent tyrosine recombinases play a major role in genome dynamics and evolution. However, despite hundreds of these systems having been identified in genome databases, very few have been described in detail, with none from phages that infect Bacillota (formerly Firmicutes). In this study, we reanalyzed the recombination module of Lactobacillus delbrueckii subsp. bulgaricus phage mv4, previously considered atypical compared with classical systems. Our results reveal that mv4 integrase is a 369 aa protein with all the structural hallmarks of recombinases from the Tn916 family and that it cooperatively interacts with its recombination sites. Using randomized DNA libraries, NGS sequencing, and other molecular approaches, we show that the 21-bp core-attP and attB sites have structural similarities to classical systems only if considering the nucleotide degeneracy, with two 7-bp inverted regions corresponding to mv4Int core-binding sites surrounding a 7-bp strand-exchange region. We also examined the different compositional constraints in the core-binding regions, which define the sequence space of permissible recombination sites.

Protein-DNA binding sites prediction based on pre-trained protein language model and contrastive learning.

Protein-DNA interaction is critical for life activities such as replication, transcription and splicing. Identifying protein-DNA binding residues is essential for modeling their interaction and downstream studies. However, developing accurate and efficient computational methods for this task remains challenging. Improvements in this area have the potential to drive novel applications in biotechnology and drug design. In this study, we propose a novel approach called Contrastive Learning And Pre-trained Encoder (CLAPE), which combines a pre-trained protein language model and the contrastive learning method to predict DNA binding residues. We trained the CLAPE-DB model on the protein-DNA binding sites dataset and evaluated the model performance and generalization ability through various experiments. The results showed that the area under ROC curve values of the CLAPE-DB model on the two benchmark datasets reached 0.871 and 0.881, respectively, indicating superior performance compared to other existing models. CLAPE-DB showed better generalization ability and was specific to DNA-binding sites. In addition, we trained CLAPE on different protein-ligand binding sites datasets, demonstrating that CLAPE is a general framework for binding sites prediction. To facilitate the scientific community, the benchmark datasets and codes are freely available at https://github.com/YAndrewL/clape.

The bacteriophage lambda integrase catalytic domain can be modified to act with the regulatory domain as a recombination-competent binary recombinase.

Site-specific recombinase Int mediates integration of the bacteriophage λ genome into the Escherichia coli chromosome. Integration occurs once the Int tetramer, assisted by the integration host factor IHF, forms the intasome, a higher order structure, within which Int, a heterobivalent protein, interacts with two nonhomologous DNA sequences: the core recombination sites and the accessory arm sites. The binding to these sites is mediated by the catalytic C-terminal domain (CTD) and the regulatory N-terminal domain (NTD) of Int, respectively. Within Int, the NTD can activate or inhibit the recombination activity of the CTD depending on whether the NTD is bound to the arm sites. The CTD alone cannot mediate recombination, and even when the NTD and the CTD are mixed together as individual polypeptides, the NTD cannot trigger recombination in the CTD. In this work, we set to determine what modifications can unlock the recombination activity in the CTD alone and how the CTD can be modified to respond to recombination-triggering signals from the NTD. For this, we performed a series of genetic analyses, which showed that a single mutation that stabilizes the CTD on DNA, E174K, allows the CTD to recombine the core DNA sequences. When the NTD is paired with the CTD (E174K) that also bears a short polypeptide from the C terminus of the NTD, the resulting binary Int can recombine arm-bearing substrates. Our results provide insights into the molecular basis of the regulation of the Int activity and suggest how binary recombinases of the integrase type can be engineered.

High-speed AFM imaging reveals DNA capture and loop extrusion dynamics by cohesin-NIPBL.

3D chromatin organization plays a critical role in regulating gene expression, DNA replication, recombination, and repair. While initially discovered for its role in sister chromatid cohesion, emerging evidence suggests that the cohesin complex (SMC1, SMC3, RAD21, and SA1/SA2), facilitated by NIPBL, mediates topologically associating domains and chromatin loops through DNA loop extrusion. However, information on how conformational changes of cohesin-NIPBL drive its loading onto DNA, initiation, and growth of DNA loops is still lacking. In this study, high-speed atomic force microscopy imaging reveals that cohesin-NIPBL captures DNA through arm extension, assisted by feet (shorter protrusions), and followed by transfer of DNA to its lower compartment (SMC heads, RAD21, SA1, and NIPBL). While binding at the lower compartment, arm extension leads to the capture of a second DNA segment and the initiation of a DNA loop that is independent of ATP hydrolysis. The feet are likely contributed by the C-terminal domains of SA1 and NIPBL and can transiently bind to DNA to facilitate the loading of the cohesin complex onto DNA. Furthermore, high-speed atomic force microscopy imaging reveals distinct forward and reverse DNA loop extrusion steps by cohesin-NIPBL. These results advance our understanding of cohesin by establishing direct experimental evidence for a multistep DNA-binding mechanism mediated by dynamic protein conformational changes.

Structural and Dynamic Features of the Recognition of 8-oxoguanosine Paired with an 8-oxoG-clamp by Human 8-oxoguanine-DNA Glycosylase.

8-oxoguanine (oxoG) is formed in DNA by the action of reactive oxygen species. As a highly mutagenic and the most common oxidative DNA lesion, it is an important marker of oxidative stress. Human 8-oxoguanine-DNA glycosylase (OGG1) is responsible for its prompt removal in human cells. OGG1 is a bifunctional DNA glycosylase with N-glycosylase and AP lyase activities. Aspects of the detailed mechanism underlying the recognition of 8-oxoguanine among numerous intact bases and its subsequent interaction with the enzyme's active site amino acid residues are still debated. The main objective of our work was to determine the effect (structural and thermodynamic) of introducing an oxoG-clamp in model DNA substrates on the process of 8-oxoG excision by OGG1. Towards that end, we used DNA duplexes modeling OGG1-specific lesions: 8-oxoguanine or an apurinic/apyrimidinic site with either cytidine or the oxoG-clamp in the complementary strand opposite to the lesion. It was revealed that there was neither hydrolysis of the N-glycosidic bond at oxoG nor cleavage of the sugar-phosphate backbone during the reaction between OGG1 and oxoG-clamp-containing duplexes. Possible structural reasons for the absence of OGG1 enzymatic activity were studied via the stopped-flow kinetic approach and molecular dynamics simulations. The base opposite the damage was found to have a critical effect on the formation of the enzyme-substrate complex and the initiation of DNA cleavage. The oxoG-clamp residue prevented the eversion of the oxoG base into the OGG1 active site pocket and impeded the correct convergence of the apurinic/apyrimidinic site of DNA and the attacking nucleophilic group of the enzyme. An obtained three-dimensional model of the OGG1 complex with DNA containing the oxoG-clamp, together with kinetic data, allowed us to clarify the role of the contact of amino acid residues with DNA in the formation of (and rearrangements in) the enzyme-substrate complex.

Analyzing paramagnetic NMR data on target DNA search by proteins using a discrete-state kinetic model for translocation.

Before reaching their targets, sequence-specific DNA-binding proteins nonspecifically bind to DNA through electrostatic interactions and stochastically change their locations on DNA. Investigations into the dynamics of DNA-scanning by proteins are nontrivial due to the simultaneous presence of multiple translocation mechanisms and many sites for the protein to nonspecifically bind to DNA. Nuclear magnetic resonance (NMR) spectroscopy can provide information about the target DNA search processes at an atomic level. Paramagnetic relaxation enhancement (PRE) is particularly useful to study how the proteins scan DNA in the search process. Previously, relatively simple two-state or three-state exchange models were used to explain PRE data reflecting the target search process. In this work, using more realistic discrete-state stochastic kinetics models embedded into an NMR master equation, we analyzed the PRE data for the HoxD9 homeodomain interacting with DNA. The kinetic models that incorporate sliding, dissociation, association, and intersegment transfer can reproduce the PRE profiles observed at some different ionic strengths. The analysis confirms the previous interpretation of the PRE data and shows that the protein's probability distribution among nonspecific sites is nonuniform during the target DNA search process.

Harnessing sulfur-binding domains to separate Sp and Rp isomers of phosphorothioate oligonucleotides.

Chemical synthesis of phosphoromonothioate oligonucleotides (PS-ONs) is not stereo-specific and produces a mixture of Rp and Sp diastereomers, whose disparate reactivity can complicate applications. Although the current methods to separate these diastereomers which rely on chromatography are constantly improving, many Rp and Sp diastereomers are still co-eluted. Here, based on sulfur-binding domains that specifically recognize phosphorothioated DNA and RNA in Rp configuration, we developed a universal separation system for phosphorothioate oligonucleotide isomers using immobilized SBD (SPOIS). With the scalable SPOIS, His-tagged SBD is immobilized onto Ni-nitrilotriacetic acid-coated magnetic beads to form a beads/SBD complex, Rp isomers of the mixture can be completely bound by SBD and separated from Sp isomers unbound in liquid phase, then recovered through suitable elution approach. Using the phosphoromonothioate single-stranded DNA as a model, SPOIS separated PS-ON diastereomers of 4 nt to 50 nt in length at yields of 60-90% of the starting Rp isomers, with PS linkage not locating at 5' or 3' end. Within this length range, PS-ON diastereomers that co-eluted in HPLC could be separated by SPOIS at yields of 84% and 89% for Rp and Sp stereoisomers, respectively. Furthermore, as each Rp phosphorothioate linkage can be bound by SBD, SPOIS allowed the separation of stereoisomers with multiple uniform Sp configurations for multiple phosphorothioate modifications. A second generation of SPOIS was developed using the thermolabile and non-sequence-specific SBDPed, enabling fast and high-yield recovery of PS substrate stereoisomers for the DNAzyme Cd16 and further demonstrating the efficiency of this method. KEY POINTS: • SPOIS allows isomer separations of the Rp and Sp isomers co-eluted on HPLC. • SPOIS can obtain Sp isomers with 5 min and Rp in 20 min from PS-ON diastereomers. • SPOIS was successfully applied to separate isomers of PS substrates of DNAzyme.

Structure and Functions of HMGB3 Protein.

HMGB3 protein belongs to the group of HMGB proteins from the superfamily of nuclear proteins with high electrophoretic mobility. HMGB proteins play an active part in almost all cellular processes associated with DNA-repair, replication, recombination, and transcription-and, additionally, can act as cytokines during infectious processes, inflammatory responses, and injuries. Although the structure and functions of HMGB1 and HMGB2 proteins have been intensively studied for decades, very little attention has been paid to HMGB3 until recently. In this review, we summarize the currently available data on the molecular structure, post-translational modifications, and biological functions of HMGB3, as well as the possible role of the ubiquitin-proteasome system-dependent HMGB3 degradation in tumor development.

Pleiotropic Effects of PhaR Regulator in Bradyrhizobium diazoefficiens Microaerobic Metabolism.

Bradyrhizobium diazoefficiens can live inside soybean root nodules and in free-living conditions. In both states, when oxygen levels decrease, cells adjust their protein pools by gene transcription modulation. PhaR is a transcription factor involved in polyhydroxyalkanoate (PHA) metabolism but also plays a role in the microaerobic network of this bacterium. To deeply uncover the function of PhaR, we applied a multipronged approach, including the expression profile of a phaR mutant at the transcriptional and protein levels under microaerobic conditions, and the identification of direct targets and of proteins associated with PHA granules. Our results confirmed a pleiotropic function of PhaR, affecting several phenotypes, in addition to PHA cycle control. These include growth deficiency, regulation of carbon and nitrogen allocation, and bacterial motility. Interestingly, PhaR may also modulate the microoxic-responsive regulatory network by activating the expression of fixK2 and repressing nifA, both encoding two transcription factors relevant for microaerobic regulation. At the molecular level, two PhaR-binding motifs were predicted and direct control mediated by PhaR determined by protein-interaction assays revealed seven new direct targets for PhaR. Finally, among the proteins associated with PHA granules, we found PhaR, phasins, and other proteins, confirming a dual function of PhaR in microoxia.

DNA methyltransferase DNMT3A forms interaction networks with the CpG site and flanking sequence elements for efficient methylation.

Specific DNA methylation at CpG and non-CpG sites is essential for chromatin regulation. The DNA methyltransferase DNMT3A interacts with target sites surrounded by variable DNA sequences with its TRD and RD loops, but the functional necessity of these interactions is unclear. We investigated CpG and non-CpG methylation in a randomized sequence context using WT DNMT3A and several DNMT3A variants containing mutations at DNA-interacting residues. Our data revealed that the flanking sequence of target sites between the -2 and up to the +8 position modulates methylation rates >100-fold. Non-CpG methylation flanking preferences were even stronger and favor C(+1). R836 and N838 in concert mediate recognition of the CpG guanine. R836 changes its conformation in a flanking sequence-dependent manner and either contacts the CpG guanine or the +1/+2 flank, thereby coupling the interaction with both sequence elements. R836 suppresses activity at CNT sites but supports methylation of CAC substrates, the preferred target for non-CpG methylation of DNMT3A in cells. N838 helps to balance this effect and prevent the preference for C(+1) from becoming too strong. Surprisingly, we found L883 reduces DNMT3A activity despite being highly conserved in evolution. However, mutations at L883 disrupt the DNMT3A-specific DNA interactions of the RD loop, leading to altered flanking sequence preferences. Similar effects occur after the R882H mutation in cancer cells. Our data reveal that DNMT3A forms flexible and interdependent interaction networks with the CpG guanine and flanking residues that ensure recognition of the CpG and efficient methylation of the cytosine in contexts of variable flanking sequences.

Residues located in the primase domain of the bacteriophage T7 primase-helicase are essential for loading the hexameric complex onto DNA.

The T7 primase-helicase plays a pivotal role in the replication of T7 DNA. Using affinity isolation of peptide-nucleic acid crosslinks and mass spectrometry, we identify protein regions in the primase-helicase and T7 DNA polymerase that form contacts with the RNA primer and DNA template. The contacts between nucleic acids and the primase domain of the primase-helicase are centered in the RNA polymerase subdomain of the primase domain, in a cleft between the N-terminal subdomain and the topoisomerase-primase fold. We demonstrate that residues along a beta sheet in the N-terminal subdomain that contacts the RNA primer are essential for phage growth and primase activity in vitro. Surprisingly, we found mutations in the primase domain that had a dramatic effect on the helicase. Substitution of a residue conserved in other DnaG-like enzymes, R84A, abrogates both primase and helicase enzymatic activities of the T7 primase-helicase. Alterations in this residue also decrease binding of the primase-helicase to ssDNA. However, mass photometry measurements show that these mutations do not interfere with the ability of the protein to form the active hexamer.

An ELISA-based platform for rapid identification of structure-dependent nucleic acid-protein interactions detects novel DNA triplex interactors.

Unusual nucleic acid structures play vital roles as intermediates in many cellular processes and, in the case of peptide nucleic acid (PNA)-mediated triplexes, are leveraged as tools for therapeutic gene editing. However, due to their transient nature, an understanding of the factors that interact with and process dynamic nucleic acid structures remains limited. Here, we developed snapELISA (structure-specific nucleic acid-binding protein ELISA), a rapid high-throughput platform to interrogate and compare up to 2688 parallel nucleic acid structure-protein interactions in vitro. We applied this system to both triplex-forming oligonucleotide-induced DNA triplexes and DNA-bound PNA heterotriplexes to describe the identification of previously known and novel interactors for both structures. For PNA heterotriplex recognition analyses, snapELISA identified factors implicated in nucleotide excision repair (XPA, XPC), single-strand annealing repair (RAD52), and recombination intermediate structure binding (TOP3A, BLM, MUS81). We went on to validate selected factor localization to genome-targeted PNA structures within clinically relevant loci in human cells. Surprisingly, these results demonstrated XRCC5 localization to PNA triplex-forming sites in the genome, suggesting the presence of a double-strand break intermediate. These results describe a powerful comparative approach for identifying structure-specific nucleic acid interactions and expand our understanding of the mechanisms of triplex structure recognition and repair.

Modeling the zinc effect on OsASR5-STAR1 promoter interaction by molecular dynamics.

Intrinsically disordered proteins (IDPs) have numerous dynamic conformations. Given the difficulties in tracking temporarily folded states of this kind of protein, methods such as molecular modeling and molecular dynamics (MD) simulations make the process less costly, less laborious, and more detailed. Few plant IDPs have been characterized so far, such as proteins from the Abscisic acid, Stress and Ripening (ASR) family. The present work applied, for the first time, the two above-mentioned tools to test the feasibility of determining a three-dimensional transition model of OsASR5 and to investigate the relationship between OsASR5 and zinc. We found that one of OsASR5's conformers contains α-helices, turns, and loops and that the metal binding resulted in a predominance of α-helix. This stability is possibly imperative for the transcription factor activity. The promoter region of a sugar transporter was chosen to test this hypothesis and free energy calculations showed how the ion is mandatory for this complex formation. The results produced here aim to clarify which conformation the protein in the bound state assumes and which residues are involved in the process, besides developing the understanding of how the flexibility of these proteins can contribute to the response to environmental stresses.

Understanding the Target Search by Multiple Transcription Factors on Nucleosomal DNA.

The association of multiple Transcription Factors (TFs) in the cis-regulatory region is imperative for developmental changes in eukaryotes. The underlying process is exceedingly complex, and it is not at all clear what orchestrates the overall search process by multiple TFs. In this study, by developing a theoretical model based on a discrete-state stochastic approach, we investigated the target search mechanism of multiple TFs on nucleosomal DNA. Experimental kinetic rate constants of different TFs are taken as input to estimate the Mean-First-Passage time to recognize the binding motifs by two TFs on a dynamic nucleosome model. The theory systematically analyzes when the TFs search their binding motifs hierarchically and when simultaneously by proceeding via the formation of a protein-protein complex. Our results, validated by extensive Monte Carlo simulations, elucidate the molecular basis of the complex target search phenomenon of multiple TFs on nucleosomal DNA.