Insights into protein–DNA interactions from hydrogen bond energy‐based comparative protein–ligand analyses

Abstract Hydrogen bonds play important roles in protein folding and protein–ligand interactions, particularly in specific protein–DNA recognition. However, the distributions of hydrogen bonds, especially hydrogen bond energy (HBE) in different types of protein–ligand complexes, is unknown. Here we performed a comparative analysis of hydrogen bonds among three non‐redundant datasets of protein–protein, protein–peptide, and protein–DNA complexes. Besides comparing the number of hydrogen bonds in terms of types and locations, we investigated the distributions of HBE. Our results indicate that while there is no significant difference of hydrogen bonds within protein chains among the three types of complexes, interfacial hydrogen bonds are significantly more prevalent in protein–DNA complexes. More importantly, the interfacial hydrogen bonds in protein–DNA complexes displayed a unique energy distribution of strong and weak hydrogen bonds whereas majority of the interfacial hydrogen bonds in protein–protein and protein–peptide complexes are of predominantly high strength with low energy. Moreover, there is a significant difference in the energy distributions of minor groove hydrogen bonds between protein–DNA complexes with different binding specificity. Highly specific protein–DNA complexes contain more strong hydrogen bonds in the minor groove than multi‐specific complexes, suggesting important role of minor groove in specific protein–DNA recognition. These results can help better understand protein–DNA interactions and have important implications in improving quality assessments of protein–DNA complex models.

involved in mediating signal transduction, regulating apoptotic pathways, and immune responses. [4][5][6] Protein-protein (PP) interactions form essential complexes like hormone-receptor, antibody-antigen, and protease-inhibitor, which control cell signaling, electron transport, signal transduction, and cell metabolism. 7 Disruptions in these interactions can cause serious medical conditions such as cancer, cardiovascular, and neurodegenerative disorders. [7][8][9][10] Knowledge of detailed interactions among these complexes at atomic resolution is therefore essential to understanding the underlying mechanisms that govern biochemical processes. It also has important implications in biomedical applications such as protein-ligand docking, in silico design of inhibitors and interfaces, and virtual screening of drugs library in the pharmaceutical industry.
Hydrogen bonds (HBs) play key roles in conferring binding specificity of macromolecular complexes. [11][12][13][14] An HB is generally considered as a weak, electrostatic interaction between a polar acceptor atom that carries a lone pair of electrons and a hydrogen atom that is covalently linked to a polar atom, oriented toward each other at an equilibrium distance. This orientation-and distance-dependent nature of HBs is vital in providing the shape and chemical complementarity for selective recognition and binding of complexes. 12 In PD complexes, for example, HBs play a key role in DNA base readout by proteins and act as the major contributor to binding specificity that is vital for the biomolecular function of protein-DNA complexes. 15 The recognition of DNA by proteins is guided by an innate hydrogenbonding pattern that generates an initial unstable nonspecific, intermediate complex with high energy. [16][17][18][19] While most of this recognition is expected to occur through the signature hydrogen-bonding pattern in the major groove, many DNA binding proteins also bind to the minor groove through hydrogen bonding and shape readout. 15,20 Later, this complex transitions to a stable and highly specific low energy state through reversible structural deformations that are also guided by a specific HB pattern. 12 In PP complexes, HBs influence stability as well as binding specificity at the interface. 14 Interfacial hydrophilic sidechains of a PP complex have a high charge density that is stabilized primarily through hydrogen bonding. Buried polar atoms at the interface not involved in hydrogen bonding may destabilize the complex. [21][22][23][24] Peptide binding, on the other hand, utilizes HBs to improve interface packing density as well as minimize the entropic cost of transitioning from a highly flexible, unstructured peptide to a well-defined rigid structure in a complex with protein. 25 On average, PT interface contains more HBs per 100 Å 2 interface area when compared to PP interface and PT interface HBs generally are more linearly oriented. 25 In addition to binding, HBs are the primary driving force in folding of protein chains into core secondary structures such as alpha helices and beta sheets and base pairing in nucleic acids. 11 HBs also bring flexibility to the structure, which is central to the dynamic nature of proteins and plays a key role in allosteric, catalytic, and binding activities. 11,26 The role of HBs in binding and folding of complexes has previously been studied as individual cases as well as a group of cases. 18,[27][28][29][30] Mandel-Gutfreund et al. studied different types of HBs at the interface of 28 X-ray crystal structures of protein-DNA complexes. The HBs were classified according to the types of donor and acceptor atoms, such as backbone, sidechain, or base edges. 13 Xu et al. performed a similar analysis on 319 protein-protein complexes. 14 London et al. compared the types of HBs at the interface and within protein chains of 103 protein-peptide complexes. They further compared the types of HBs in protein-peptide complexes to those in protein-protein complexes. 25 Rawat and Biswas in 2011 performed a comparison of HBs along with several other structural features to investigate the role of flexibility in protein-DNA, protein-RNA, and protein-protein complexes. 31 Jiang et al. demonstrated that in protein-protein complexes, the average energy contribution of a HB is $30%. 32 Zhou and Wang recently compared short HBs, where donor-acceptor distance is less than 2.7 Å, in 1663 high-quality protein, protein-ligand, and protein-nucleic acid structures. 33 Itoh et al. showed that the interaction energy of even the weaker N + C HÁÁÁO HBs is comparable to other protein-ligand interactions such as π-π interactions suggesting the importance of considering HB energy in drug design. 34 While analyses based on the number of HBs with a single energy cutoff or a distance/angle cutoff can provide useful information about the role of HBs in protein-ligand interaction, they have an intrinsic flaw since strong and weak HBs are treated equally. Moreover, the distributions of interfacial HBs in terms of HB strength or HB energy in protein-ligand complexes, and more importantly, the distributions of interfacial HB energy among different types of protein-ligand complexes remain unknown. To address these issues, in this study, we performed a holistic statistical comparative analysis of HBs across interfaces and within protein chains (intrachain) among PP, PT, and PD complexes to get an insight into their roles in each type of complexes. In addition to comparing the types and locations of HBs in each type of complexes, we investigated the HB energy distributions and found significant differences among these three types of complexes, especially a unique pattern in protein-DNA complexes. To the best of our knowledge, an HB energy-based large-scale comparison of macromolecular complexes has never been explored before.
Since the M-TASSER dimer library was published over 10 years ago, we generated an updated dataset, called Protein Homo/ Heterodimer Library (PHDL) using some of the guidelines described in the original paper (Table S1).
Each of the three datasets for PD represents a specific category of protein-DNA complexes. The HS dataset comprises 29 PD complexes with high binding specificity between protein and DNA whereas the MS dataset comprises 104 cases, in which proteins can bind to multiple conserved DNA sequences. 3 The RDPD dataset consists of 38 highly diverse nonredundant TF-DNA complexes that cover 11 structural folds, 15 super-families, and 28 families. 35 The two PT complex datasets differ mainly in the peptide chain lengths. InterPep comprises protein complexes with peptides ranging from 5 to 25 amino acids whereas peptides in LEADS-PEP are 3-12 amino acids long. 36,37 InterPep is a larger dataset with 502 X-ray and NMR structures, which was originally developed for testing a peptidebinding site prediction pipeline. 37 LEADS-PEP, on the other hand, is a much smaller dataset with 53 carefully curated and widely used complexes designed specifically for peptide-based therapeutics and peptide docking. It contains only X-ray crystal structures with a resolution better than 2 Å. 36 The complexes in the PP datasets differ mainly in size and definition of interaction unit. The protein-protein docking benchmark (RDPP) has 230 complex structures that were experimentally solved with corresponding unbound components available. 39 The structures in the RDPP dataset represent a diverse combination of antigen-antibody, enzyme-substrate, enzyme-regulatory complex, GPCR proteins, and several other classes of proteins. The docking benchmark defines a true interaction as one that has functional significance as identified in the literature and agreed upon by the scientific community. The second PP dataset PHDL, a protein homo/heterodimer library, determines the oligomeric state from PDB files. 40 PHDL contains nonredundant heterodimers (Table S1A) and homodimers (Table S1B), where no two chains share more than 30% sequence identity with each other and each interacting partner has at least 40 amino acids.
In addition to these individual datasets, we pooled the datasets of the same type of complexes together and generated three larger, nonredundant, and highly diverse datasets (

| Dataset processing
The datasets were filtered rigorously for accurate analysis. In case of multiple models for one native structure as in the NMR entries, only the first model was selected. All the heteroatoms, including water molecules, were removed since we do not consider solvation effects for the sake of simplicity and fair comparison. Proteins that have residues with insertion codes were renumbered accordingly. Since considering the alternate locations of a residue in an experimentally solved crystal structure may result in over counting the number of HBs, only the state with the highest occupancy for a given residue was included for analysis. The complexes with internal missing residues, that is, residues that are not on the N or C terminal of the chain were discarded.

| Identification of HBs
Two widely used HB annotation programs, FIRST (Floppy Inclusion and Rigid Substructure Topography) and HBPLUS, were used to identify HBs with default parameters. 42,43 Reduce was used to add hydrogen atoms to pdb files for FIRST HB calculations while HBPLUS calculates the hydrogen atom positions within the program. 44 FIRST employs an energy-based approach and the HB energy is calculated as in Equation (1). 42,45 where d is the donor-acceptor distance. d 0 (2.8 Å) and V 0 (8 kcal/mol) represent the equilibrium distance and well-depth, respectively. 46 The angle term F(θ,ϕ,φ) is calculated based on the hybridization state of the acceptor and donor atoms, where θ is the donor-hydrogenacceptor angle, ϕ is the hydrogen-acceptor-base angle, and φ is the angle between the normals of the planes defined by the six atoms attached to the sp 2 center as described by Dahiyat et al. 45 The FIRST program was used for both the number of HBs annotations using a widely used HB energy cutoff of À0.6 kcal/mol as well as for HB energy-based analysis. HBPLUS identifies HB with a distance-angle approach and defines the optimal distance between the donor and acceptor as 2.5 Å or smaller and the optimal angle as 90 or higher. 43

| Interface analysis and comparison
Since the interface sizes are different among different types of complexes (Table 1), in order to accurately assess the roles of HB at the interface of PP, PT, and PD complexes, the numbers of HBs were compared with respect to the interfacial surface area. The interfacial surface area (iSA) of a complex was calculated using NACCESS v2.1.1 with default parameters as shown in Equation (2): where SA p and SA L represent the surface area of protein and ligand, respectively, and SA C is the surface area of the protein-ligand complex. For multichain components, SA P is the surface area of the protein unit while SA L is the surface area of the ligand unit.
The HB distributions were compared at three different aspects:  Table 2.

| HBs at the interface of complexes
We first compared the number of HBall and HBSP in PDnrall, PPnrall, and PTnrall datasets. Based on HB annotations from FIRST with the widely used energy cutoff of À0.6 kcal/mol, 48 we found that the number of interface HBall and the number of interface HBSP in PD complexes are significantly higher than those in the PP and PT complexes (Figure 2A,B). The number of HBall and HBSP in PT complexes are significantly less than those in PP complexes (Figure 2A,B). Results from HBPLUS are consistent with the data from FIRST except that the number of HBSP in PP complexes is larger than that in PD complexes with HBPLUS ( Figure S1A,B). Interestingly, when the FIRST energy cutoff is set at À0.1 kcal/mol, the results are more similar to the HBPLUS data ( Figure S2A,B).
Since the interface areas among the three types of complexes are different with PP complexes having the largest average interfacial area and PT complexes having the smallest average interfacial area (Table 1), comparing the raw number of interface HBs might be biased towards the complexes with a larger contact surface. Therefore, we normalized the number of interface HBs, HBall, and HBSP, by the interfacial surface area (iSA). Figure 2C,D shows that both HBall/iSA and HBSP/iSA ratios of PD complexes are significantly higher than those in the PP complexes and PT complexes. There is a clear pattern for the iSA normalized HBSP, PD > PP > PT. When the analyses were carried out with HBPLUS, the results are consistent with the results from FIRST ( Figure S1). Even though no significant difference of the ratio HBall/iSA from FIRST is found between PP and PT complexes for a two-tailed test ( Figure 2C), one-tailed test with a null hypothesis that HBall/iSA in PP is not smaller than HBall/iSA in PT results a p- area than that in PP. 25 However, the HBSP/iSA ratio is the opposite, suggesting relatively fewer interface HBSP in PT complexes when the interface area is taken into consideration.

| Types of HBs at interface and within intrachain
We compared the distributions of the HB types at complex interface or within protein (intrachain) in PP, PT, and PD complexes and between individual complexes of the same type of complexes. Figure 3A and Table 3 show that there is no significant difference among the types of HBs within proteins in all three types of complexes. BB-BB HBs represent the largest number of overall HBs within proteins (66%-69%) followed by the Mixed (17%-20%), and SC-SC (14%) HBs, respectively ( Figure 3A). This is not surprising because the two major secondary structure types of the core protein structure, α-helices, and β-sheets, are stabilized by backbonebackbone HBs.
The distributions of the HB types at interface, however, are significantly different from the intrachain and among the three types of complexes ( Figure 3B and Table 3 Figures S5B and S6).

| Strength of HBs at interface and within protein chain
We classified the strength of HBs into four categories based on HB energy from the FIRST program with different energy cutoffs used in previous studies as shown in Table 2. 17,42,48 For intrachain HBs within proteins, no significant differences were found among the three types of complexes ( Figure 4A and However, the interface HB energy distributions among different types of complexes are significantly different and exhibit a unique pattern for the PD complexes ( Figure 4B and Table 4). There  (Table 4).
We also compared the energy distributions of each HB type across interfaces ( Figure 5B). Similar to the pattern observed for all HBs in PD, energy distributions of different types of interface HB in PD complexes also differ significantly from PP and PT complexes while there is no significant difference between PP and PT complexes (Table S3)

| Comparison of HBs between HS and MS datasets
In our previous study, we demonstrated that highly specific HS protein-DNA complexes have more HBs than the multi-specific MS protein-DNA complexes, including both total HBs and sidechain-base HBs. 3 It is intriguing to see whether there is any relationship between the HB strength and protein-DNA binding specificity. We first compared the HB types and energy categories within proteins as well as at the interface of HS and MS complexes. No significant differences between HS and MS complexes were found in terms of energy categories ( Figures S7 and S8) while there are significant differences between the intrachain and interface for both HS (p-value: 9.673eÀ07) and MS complexes (p-value: 6.413eÀ07). We did observe some statistically nonsignificant small differences. For example, the number of SC-SC interface HBs in HS (32%) is slightly higher than that in MS (28.2%; Figure S5A). Both HS and MS complexes show similar interface HB energy distributions with an overall balance of strong and weak HBs, but HS complexes have a slightly higher percentage of HBs in category IV ( Figure S7).
Since both major and minor grooves are known to play important roles in the base and shape readout mechanisms in specific protein-DNA recognition, 3,15,20,49 we compared the energy distributions of total HBs and sidechain-base HBs in the major and minor grooves. Our comparative analyses on energy categories are based on HB energy cutoffs (À0.1, À0.6, À1.0, and À1.5 kcal/mol) from previous studies ( Table 2). 17,42,48 To test if similar results can be observed with different HB energy discretization, the HBs were grouped using a larger energy range separated by À0.1, À0.7, À1.3, and À2.0 kcal/mol (  (Figures 4-7, Tables 4 and S3) with energy ranges in Table 2, suggesting our key findings are not affected by different discretization of HB energy.
The above findings have important functional and practical implications. While omitting HB information in assessing predicted PP and PT complex models may have minimal effect, our results suggest consideration of HBs is beneficial to quality assessment of protein-DNA complex models since both the raw number and the normalized number of interface HBs in PD complexes are much higher than those in PP and PT complexes. The use of conserved numbers of native HBs in models was suggested to evaluate the quality of protein-peptide models. 50 We found that using the number of HBs can improve quality assessment of protein-DNA complex models. 51 However, due to the unique pattern of interface HB energy distributions in PD complexes and the dynamic nature of macromolecules, it could help model evaluations by considering the HB energy instead of using the raw number of HBs. We demonstrated in our previous study that the accuracy of structure-based prediction of transcription factor binding sites could be improved by adding an HB energy term. 52,53 Our data also provide an insight into the mechanism of binding specificity between protein and DNA. We observed an approximate balance of high and low energy interface HBs in PD complexes, but not in the other two types of complexes ( Figures 4B and 5B). One possibility of such difference lies in the geometry of interacting components as geometry is one of the key factors affecting HBE and strength. 46 While DNA is not a rigid molecule, the double-helical nature restricts the atoms that can form optimal HBs with protein sidechains while the peptide and protein surfaces have a relatively higher flexibility to position atoms for stronger HBs. Other than the unique structure of DNA double-helix that contributes to the pattern of energy distribution, it may also reflect the kinetics of protein-DNA recognition and binding, and the functions of many DNA binding pro- Highly specific DNA binding proteins have more strong HBs than the MS group comprising transcription factors (Figures 6 and 7). 3 The most interesting finding is from the DNA minor groove HB analysis. Both the energy of all HBs and the sidechain-base HBs of highly specific protein-DNA complexes are significantly different from that of multi-specific protein-DNA complexes (Figure 7). While it is generally thought that minor groove contacts play little role in conferring specific protein-DNA interactions, more studies have shown that this might not be the case. It has been reported that local sequence-dependent minor groove shape plays an important role in specific recognition between protein and DNA. 15,20,[55][56][57] The number of contacts in minor grooves of HS complexes is more than that in MS complexes and the HS complexes contain wider minor grooves than MS, 3 thus making it possible for optimal orientation of atoms to form stronger HBs. Our results further demonstrate that the minor groove HBs play more critical roles in conferring binding specificity than previously thought.

CONFLICT OF INTEREST
The authors have declared no conflicts of interest for this article.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available in the supplementary file and references cited in the paper. The PDB files are from https://www.rcsb.org.