How Divalent Cations Interact with the Internal Channel Site of Guanine Quadruplexes.

The formation of guanine quadruplexes (GQ) in DNA is crucial in telomere homeostasis and regulation of gene expression. Pollution metals can interfere with these DNA superstructures upon coordination. In this work, we study the affinity of the internal GQ channel site towards alkaline earth metal (Mg2+ , Ca2+ , Sr2+ , and Ba2+ ), and (post-)transition metal (Zn2+ , Cd2+ , Hg2+ , and Pb2+ ) cations using density functional theory computations. We find that divalent cations generally bind to the GQ cavity with a higher affinity than conventional monovalent cations (e. g. K+ ). Importantly, we establish the nature of the cation-GQ interaction and highlight the relationship between ionic and nuclear charge, and the electrostatic and covalent interactions. The covalent interaction strength plays an important role in the cation affinity and can be traced back to the relative stabilization of cations' unoccupied atomic orbitals. Overall, our findings contribute to a deeper understanding of how pollution metals could induce genomic instability.

Introduction

Chronic or acute exposure to pollution metals, such as mercury, cadmium, and lead, is known to have negative effects on the health of humans and other organisms as these elements can cause genomic instability. One possible mechanism relevant to the observed genotoxic effects is the coordination of pollution metal cations to guanine quadruplex DNA. These structures are crucial in terms of gene expression and genome stability as they are located in telomeric and gene promotor regions of DNA, that have high concentrations of guanine residues. If specific environmental conditions are met, these guanine‐rich sequences can fold into more complex structures, called guanine quadruplexes (GQs, Figure 1). The basic building block of a GQ consists of four guanine bases interconnected by Hoogsteen hydrogen bonding, a so‐called guanine quartet (G4).

Figure 1

Schematic representation of the examined double‐layer guanine quadruplexes (GQ) hosting a metal cation M . The structures are shown of the guanosine phosphate dimer as building block of the GQ (left), and a guanine quartet (G4) with hydrogen bond labelling (right).

In our previous work on GQs, we showed that the hydrogen bonds in G4 experience a large synergistic effect. This cooperativity in G4 originates from charge separation occurring through donor‐acceptor interactions in the σ‐electron system. Additional layers of guanine quartets can stack on top of each other, where the four lateral pillars are constituted by a sugar‐phosphate backbone (Figure 1). The several different kinds of loops interconnecting these pillars contribute to the very large topological variety of these structures. The carbonylic oxygens of the guanine bases define a partially negatively charged channel at the very core of GQs that under biological conditions hosts the monovalent alkali cations K+ or Na+ (Figure 1). It was demonstrated that the coordination of alkali metal cations to GQs is essential to overcome the entropic penalty of self‐assembly and plays a crucial role in maintaining GQ stability even in the presence of damaged nucleobases. The cation‐quadruplex interaction follows an order of affinity determined by the size of the cation and its desolvation energy, as the cations need to be desolvated completely upon entering the negative GQ channel. Nieuwland et al. showed that the alkali metal cation affinity order observed in double‐layer GQs (K+>Na+>Rb+>Li+) is preserved in multi‐layer systems as they occur in vivo.

Although it is well established that divalent metal cations (i. e. ionic charge of +2) generally bind GQ structures more strongly than the monovalent ones (i. e. ionic charge of +1), the order of their affinity has no experimental consensus nor a fully understood theoretical explanation. Even though evidence of GQ formation in the presence of divalent cations is available,[ , ] our computational results are not meant to rule out the formation of alternative structures, like metal‐mediated base pairs, DNA‐DNA adducts as triplexes, or other non‐canonical motifs, because DNA is highly versatile and divalent cations can induce strong non‐specific interactions.

In this work, we aim at getting a deeper understanding of how divalent cations could induce genome instability, by studying the affinity of the GQ internal channel site towards a wide set of cations using dispersion‐corrected density functional theory (DFT‐D) based computations with implicit aqueous solvation. Our quantum chemical computations analyze the interaction of GQs with divalent cations of alkaline earth metals (Group 2: Mg2+, Ca2+, Sr2+, and Ba2+), transition metals (Group 12: Zn2+, Cd2+, and Hg2+), and a post‐transition metal (Group 14: Pb2+). We want to highlight the choice of the Cd2+, Hg2+, and Pb2+ cations since they represent some of the most widespread heavy metal pollutants. We have quantified and analyzed the energetic parameters that determine the order of the GQ‐cation affinity using the state‐of‐the‐art activation strain model (ASM) followed by an energy decomposition analysis (EDA) and quantitative Kohn‐Sham molecular orbital analysis. We highlight the relationship between cation desolvation, the size of cations, and the strength of interaction in the observed cation affinities. Furthermore, we demonstrate how most of the divalent cations could potentially displace K+ in vivo and could lead to genotoxicity due to structural perturbation of telomeric and promotor regions.

Computational Details

All calculations were performed with the Amsterdam Density Functional (ADF, version 2016.102 for GQs, 2017.208 for the formaldehyde model systems) program using dispersion‐corrected relativistic density functional theory (DFT‐D) at the ZORA‐BLYP‐D3(BJ)/TZ2P level of theory for geometry optimizations and stationary point energies.[ , ]

Solvent effects were accounted for by using the conductor‐like screening model (COSMO, using the Esurf cavity construction), as implemented in the ADF program. The COSMO solvation model has proven reliable for the optimization of DNA in water and a spectrum of other solvents. Radii of cations have been computed according to the procedure presented in Ref. [18b]. The solvation radii parameters used to generate the COSMO cavities in ADF have been systematically tuned so that the computed electronic solvation energy of the cations matches with the experimental cation Gibbs free energy of solvation (see SI Table S1). The energies of the empty GQ scaffold and the GQ‐M complexes correspond to the electronic solvation energies computed using COSMO since computation of Gibbs free energies of biomolecular systems of such significant size is only feasible at lower, less accurate levels of theory. We chose this method as this work aims at elucidating cation binding affinity trends, rather than reproducing experimental values. For an estimate of the Gibbs free energy of formation of an empty GQ and GQ‐K+ we refer the reader to our previous work in Ref. [6a].

No geometrical or symmetry constraints have been imposed on the GQs nor guanosine dimers. The quadruplex structures in this work are all parallel‐stranded right‐handed GQs with anti‐glycosidic torsion angles at all guanines, meaning that all eight hydrogen bonds of the double layer point in the same direction. DNA‐GQs can exhibit various configurations.[ , ] In our previous work we showed how parallel and anti‐parallel conformation furnish the same trend in cation affinity for alkali metals and therefore will interact in a similar fashion with divalent cations.

As shown in Figure 1, the negatively charged sugar‐phosphate backbone of the investigated quadruplexes has been neutralized with the use of H+. As well as previously demonstrated, we show in the Supporting Information (SI) that the use of H+ counterions does not alter the cation affinity order, compared to the biologically occurring Na+ counterions (SI Table S2 and S3).

Cartesian coordinates and total bond energies of the optimized structures can be found in the Supporting Information. For the coordinates and energies of the structures without cations and comprising alkali metal cations, we direct the reader to our previous work in Ref. [8].

The energy of formation and the bonding between the metal cation and the GQ scaffold are the two primary measures used in this work to describe the phenomena of quadruplex association and complexation. The energy of formation (ΔE form) is based on the association of four individually optimized guanosine dimers with a cation into a guanine quadruplex (see Figure 2 and Equation 1:

Figure 2

Definition of the formation energy (ΔE form) of the GQ‐M complexes in aqueous solution, with the studied metal cations M and their corresponding groups in the periodic table of elements. The blue circles surrounding the structures indicate that the calculation is performed with COSMO to simulate solvation in water.

The term E(GQ‐M )aq denotes the energy of the optimized double‐layered GQ containing the metal cation. E(GG)aq is the energy of a single optimized guanosine dimer (Figure 1). The term E(M )aq expresses the energy of the cation M with charge n+, with n=1 or 2. The examined metal cations are listed in Figure 2. Every energy term of Equation 1 has been calculated using implicit solvation in water, indicated by the ‘aq′ subscript.

The bond energy (ΔE bond) is the energy change associated with the complexation of the empty GQ scaffold with a cation (Figure 3 and Equation 2). The term E(GQ)aq is the energy of the optimized empty double‐layered GQ scaffold in the solvated state.

Figure 3

Partitioning of the bond energy (ΔE bond) of the complexation of a metal cation M to the empty cavity of the GQ scaffold, yielding the final GQ‐M complex. The blue circles surrounding the structures indicate that the calculation is performed with COSMO to simulate solvation in water.

To get more insight into the observed bond energy trends, ΔE bond can be decomposed into various terms as formulated by Equation 3, according to the activation strain model of reactivity and bonding (see Figure 3).

The strain energy ΔE strain is calculated in the gas phase and is the energy required to deform the solvated state geometry of the empty GQ scaffold (E(GQ)gas) to the solvated state geometry it acquires when it binds the cation in the final complex (E(GQ‐[ ])gas) (Equation 4). So the geometries of GQ and GQ‐[ ] are obtained from optimizations with COSMO, but the stationary point energies are computed in the gas phase, indicated by the ‘gas’ subscript.

The desolvation and solvation energy can be computed as the energy difference between the solvated species and the stationary point energy of the solvated species in the gas phase (Equations 5 and 6.

Finally, the interaction energy ΔE int of the cation with the prepared (i. e. deformed) GQ scaffold in the gas phase is formulated by Equation 7. As before, this involves gas phase stationary point calculations of COSMO optimized geometries.

The interaction energy in this model is examined in the framework of the Kohn‐Sham molecular orbital theory using a quantitative energy decomposition analysis (EDA) that divides the total interaction (ΔE int) into electrostatic interaction (ΔV elstat), Pauli repulsion (ΔE Pauli), orbital interaction (ΔE oi), and dispersion (ΔE disp) components:

The ΔV elstat term corresponds to the classical electrostatic interaction between the unperturbed charge distributions of the prepared (i. e. deformed) GQ scaffold and M and is usually attractive. The Pauli repulsion ΔE Pauli comprises the destabilizing interaction as the result of overlapping occupied orbitals and is responsible for any steric repulsion. The orbital interaction ΔE oi accounts for the charge transfer (i. e. donor‐acceptor interactions between occupied orbitals on one molecular fragment and unoccupied orbitals on the other, including HOMO‐LUMO interactions) and polarization (empty‐occupied orbital mixing on one fragment due to the presence of the other fragment). In addition, dispersion interactions are included by the ΔE disp correction term.

Results and Discussion

Geometries and Energies of Formation of GQ‐M

The study of the influence of divalent cations on the structure and stability of guanine quadruplexes (GQs) begins with the analysis of the geometrical parameters and the energies of formation (ΔE form). Divalent cations are known to interact similarly with the internal channel site of quadruplexes as their monovalent counterparts. The calculated geometric values are reported for different M cations in Table 1 and the geometries are presented in Figure 4. We are interested in the relative stability of quadruplexes comprising divalent cations, compared to the conventional biologically occurring monovalent alkali cations. Therefore, the GQ‐K+ structure, as reported in our previous work, will be used as a reference as it generates the most stable GQ structure among the alkali metal cation series.

Table 1

Energies of formation (ΔE form in kcal mol−1) and geometrical parameters (in Å) of the optimized GQ‐M complexes.[a,b]

Group[b]	M n+	Ionic radius[c]	d[O⋅⋅⋅M n+][d]	R [e]	N2(H)⋅⋅⋅N7[f]	N1(H)⋅⋅⋅O6[g]	ΔE form
–	None	–	3.03	3.60	2.88	2.81	−62.5
1	K+	1.38	2.82	3.33	2.88	2.82	−115.4
2	Mg2+	0.72	2.05/4.04	3.23	2.87/2.81	2.80/2.86	−105.9
	Ca2+	1.00	2.54	2.69	2.83	2.85	−131.3
	Sr2+	1.13	2.64	2.86	2.86	2.86	−131.6
	Ba2+	1.36	2.77	3.13	2.90	2.87	−130.4
12	Zn2+	0.75	2.07/3.94	3.16	2.78/2.88	2.80	−121.2
	Cd2+	0.95	2.55	2.72	2.81	2.81	−119.0
	Hg2+	1.02	2.63	2.88	2.83	2.81	−133.4
14	Pb2+	1.18	2.70	3.01	2.87	2.84	−136.4

[a] Computed at COSMO(H2O)‐ZORA‐BLYP‐D3(BJ)/TZ2P. [b] Data of M =None, Group 1 taken from our previous work (see Ref. [8]). [c] Ionic radii taken from Ref. [21]. [d] Average distance between the carbonylic oxygen atoms and the metal cation. For the empty scaffold, the midpoint of the eight carbonylic oxygen atoms was taken. For Mg2+ and Zn2+, two values are presented, as the quartets are not equal, i. e., cations lie in the middle of one of the quartets. [e] Difference in average z‐coordinate of the upper and lower carbonylic oxygen atoms. [f] Average outer hydrogen bond distance N2(H)⋅⋅⋅N7. For Mg2+ and Zn2+, two values are presented, as the quartets are not equal, i. e., cations lie in the middle of one of the quartets. [g] Average inner hydrogen bond distance N1(H)⋅⋅⋅O6.

Figure 4

Structures of GQ‐M , where M is a monovalent (n=1) or divalent (n=2) cation, optimized at ZORA‐BLYP‐D3(BJ)/TZ2P using COSMO to simulate solvation in water. The formation energies are shown in brackets below the structures in kcal mol−1. GQ‐M1+ structures and energies are taken from Ref. [8].

M

d[O⋅⋅⋅M ][d]

–

–

The central cavity of an empty GQ scaffold contracts upon complexation with a metal cation. This is reflected by the O⋅⋅⋅M distance and the rise (R, difference in average z‐coordinate of the upper and lower carbonylic oxygen atoms) that decreases for all metal cations, compared to the empty GQ reference. This contraction is more pronounced for the divalent cations M2+ of Group 2 of the periodic table (Mg2+, Ca2+, Sr2+, and Ba2+), as the average O⋅⋅⋅M distance goes from 2.82 Å for K+ to 2.54 Å for Ca2+, and R decreases from 3.33 Å for K+ to 2.69 Å for Ca2+. Ca2+ has the smallest ionic radius of the Group 2 elements that lie in the center of the two quartet layers, leading to the largest shrinking of the central cavity (Table 1). Mg2+ is too small to interact with all eight oxygens of the guanine bases simultaneously and resides in the middle of one quartet (Figure 4), as is the case for Li+.[ , ]

For the cations of Group 12, the same trend was observed: smaller cations lead to higher deformation of the cavity. The R distance for GQ‐Cd2+ is diminished by 0.61 Å compared to the K+ reference (Table 1). Again, the ionic radius of Zn2+ is too small for efficient coordination to two quartets, so it adopts a position in the middle of one quartet. The GQ cavity shrinks only modestly for the Group 14 Pb2+ cation compared to the K+ reference.

The variation of the hydrogen bond lengths between the guanine bases is much smaller for the GQs hosting different cations (Table 1, and Figure 1 for hydrogen bond labelling). Although, a slight shortening of the hydrogen bond distances is observed with decreasing cationic radius. It is worth noting that the computed geometrical N2(H)⋅⋅⋅N7 and N1(H)⋅⋅⋅O6 distances are consistent with experimental data.

The energies of formation of the GQ‐M structures, ΔE form, show pronounced differences between monovalent and divalent cations. In Group 2, GQ‐Ca2+, GQ‐Sr2+, and GQ‐Ba2+ have very similar ΔE form of −131.3, −131.6, and −130.4 kcal mol−1, respectively. This is about 15 to 16 kcal mol−1 more stabilizing than for GQ‐K+. Among the Group 2 cations, Mg2+ is not likely to foster the formation of GQs with an ΔE form of only −105.9 kcal mol−1, similarly as observed for Li+ within the alkali metal series.

In Group 12, the GQ‐Zn2+ and GQ‐Cd2+ species lead to formation energies of −121.2 and −119.0 kcal mol−1, respectively. Although these ΔE form values are considerably less favorable than those of the Group 2 cations, they are still more stabilizing than the ΔE form observed for GQ‐K+. Zn2+, despite having a very similar ionic size as Mg2+ and Li+, and preferentially coordinating in‐plane with one of the quartets, has a very favorable formation energy and seems to be able to form stable GQ complexes. However, this fundamental theoretical observation might never find experimental confirmation, since other processes, like the bridging interaction of Zn2+ with N7 and O6 of the guanines, could inhibit the process of GQ formation. GQ‐Hg2+and GQ‐Pb2+ have the most stabilizing formation energies of −133.4 and −136.4 kcal mol−1, respectively. The high stability of GQs complexed with these heavy metal cations might be an explanation for the genotoxicity following human exposure to lead or mercury.

Partitioning of the Bond Energy

Partitioning the M bond energy ΔE bond, as illustrated in Figure 3, allows for a better understanding of the coordination affinity for the different metal cations to the internal channel of the GQ. The results of this analysis are presented in Figure 5 (see SI Table S2 for the numerical values). In this partitioning, solvation effects are covered by the ΔE desolv and ΔE solv terms. Cations need to be completely desolvated to enter the channel of the GQ, leading to an overall destabilizing solvent effect reflected by the sum of these two terms, ΔE desolv+solv. The effect of the size of the cation is reflected by the ΔE strain term and is related to the deformation of the GQ scaffold to facilitate cation coordination. Finally, any stabilizing interactions between the cation and the GQ channel are covered by the ΔE int term.

Figure 5

Partitioning of the bond energy (ΔE bond) of metal cation M coordination to the internal channel of the GQ scaffold generating a GQ‐M complex, calculated at ZORA‐BLYP‐D3(BJ)/TZ2P using COSMO to simulate solvation in water. The energies for GQ‐M1+ are taken from Ref. [8]. The dotted lines between the data points are there to guide the reader.

Our previous work on the GQ affinity for monovalent alkali metal cations[ , ] has shown that ΔE bond is the result of a large destabilizing desolvation effect (ΔE desolv+solv), counteracted by a large stabilizing interaction energy (ΔE int). Small cations are engaged in more favorable interactions with the GQ cavity but deal simultaneously with a higher energetic cost by desolvation of the cation (Figure 5). This means that if the interactions of the ion with the GQ cavity are strong, the desolvation is also strong. The competition between ΔE int and the desolvation of the cations leads to an optimum cation affinity (ΔE bond) for K+, but ΔE bond is only 0.9 kcal mol−1 less stabilizing for Na+. Furthermore, we established that the strain energy (ΔE strain) varies less among the different alkali cations. However, a larger ΔE strain was observed for smaller cations. This larger deformation is related to the shrinking of the GQ cavity for ions with a small radius (vide supra).

The partitioning of ΔE bond for the divalent cations M2+ is also presented in Figure 5 (see SI Table S2 for numerical values). Compared to the monovalent alkali cations, the divalent cations bind with higher affinities (i. e. more stabilizing ΔE bond) to the central GQ channel. With ΔE strain and ΔE desolv+solv becoming more destabilizing for the divalent cations, this enhanced affinity originates from more stabilizing interactions (ΔE int). What causes the more favorable ΔE int will be addressed in our energy decomposition analysis in the next section.

Among the alkaline earth metal cations (Group 2, Mg2+, Ca2+, Sr2+, and Ba2+) the bond energy decomposition terms show similar trends as for the alkali metals. The smaller Group 2 ions require in general more deformation of the GQ scaffold (ΔE strain) than the larger ones but are simultaneously involved in more stabilizing interactions (ΔE int) and have higher desolvation energies (ΔE desolv+solv). Overall, this leads to a minimum ΔE bond for Sr2+ with −69.3 kcal mol−1, followed closely by Ca2+ (−69.1 kcal mol−1).

The bond energies (ΔE bond) for the divalent (post‐)transition metal cations (Group 12: Zn2+, Cd2+ and Hg2+ and Group 14: Pb2) are also obtained from a large attractive interaction energy and large desolvation energy of the cation. We observe remarkably high affinities for Hg2+ and Pb2+ coordination, with ΔE bond values of −71.1 and −74.2 kcal mol−1, respectively. This highlights the effect of heavy metal cations on the stability of genetic structures. For the cations of Group 12 and 14, the deformation of the GQ scaffold is also larger for the smaller cations as they shrink the cavity of the quadruplex to interact with the carbonyl groups of the guanine bases. Noteworthy is the more stabilizing ΔE int and more destabilizing ΔE desolv+solv for Hg2+ compared to Cd2+. The reason for the interruption of the trend for the desolvation and interaction energy terms for transition metals compared to the alkali and alkaline earth metals will be explained in more detail by our molecular orbital analysis (vide infra).

Energy Decomposition Analysis

In this section, we aim to get a deeper understanding of the trends in the interaction energies (ΔE int) of the cations M with the central, prepared, GQ cavity. In addition, we want to trace the origin of the more favorable ΔE int for the divalent cations, compared to the monovalent ones. For this purpose, ΔE int was decomposed into physically meaningful terms by performing an energy decomposition analysis (EDA, see Equation 8). The results of this analysis are presented in Figure 6 and reported in SI Table S4.

Figure 6

Energy decomposition analysis of the interaction energy (ΔE int, in kcal mol−1) between the metal cations M and the GQ scaffold in GQ‐M , calculated at ZORA‐BLYP‐D3(BJ)/TZ2P. The energies for GQ‐M1+ are taken from Ref. [8]. The dotted lines between the data points are there to guide the reader.

From Figure 6 it becomes evident that the more favorable interaction energy (ΔE int) for the divalent metal cations compared to the alkali metal cations originates from both more stabilizing electrostatic (ΔV elstat) and orbital interactions (ΔE oi). The Pauli repulsion (ΔE Pauli) becomes more destabilizing for the divalent cations, while the dispersion energy (ΔE disp) changes only minorly. The enhanced ΔV elstat can be rationalized by the increased ionic charge from +1 to +2. We explain the origin of the more stabilizing ΔE oi in the molecular orbital analysis in the next section.

Like previously found for the monovalent alkali metals,[ , ] for the group of divalent alkaline earth metals, the strength of ΔE int is inversely proportional to the ionic radius. In the series Mg2+ to Ca2+ to Sr2+ to Ba2+, the steep weakening of ΔE int is mainly determined by less favorable orbital interactions (ΔE oi), alongside the growth of steric repulsive forces (ΔE Pauli). The electrostatic interaction (ΔV elstat) becomes slightly more stabilizing from Mg2+ to Ca2+ due to switching from in‐plane coordination to inter stack coordination but then decreases for Sr2+ and Ba2+. This latter effect can be rationalized by the increasing O⋅⋅⋅M2+ distance going from Ca2+ to Sr2+ to Ba2+ (Table 1), where a longer distance leads to a less stabilizing electrostatic interaction between the cation and the guanine bases. The dispersion energy (ΔE disp) becomes only modestly more stabilizing for the heavier elements of this group.

For the (post‐)transition metal cations, similar observations hold (Figure 6). ΔV elstat becomes less stabilizing along the trend Zn2+>Cd2>Hg2+>Pb2+, which follows the trend of the increasing O⋅⋅⋅M2+ distance. In addition, the dispersion energy component becomes more favorable for heavier elements. Due to the introduction of a full d shell, the ΔE Pauli term is relatively high for the in‐plane coordinated Zn2+ cation compared to the alkali and alkaline earth metal cations Li+ and Mg2+. After zinc, ΔE Pauli becomes less destabilizing for Cd2+, but then increases again for Hg2+ and is the most destabilizing for Pb2+. However, none of these terms dictate the trend in ΔE int, as this is determined by the trend in orbital interactions. Both the ΔE int and ΔE oi terms become less stabilizing from Zn2+ to Cd2+ but become more stabilizing for the heavier cation Hg2+, and then decrease again for Pb2+. This remarkable more stabilizing ΔE int due to more stabilizing ΔE oi for Hg2+ compared to Cd2+ will be explained in the next section.

From our bond energy analysis, particularly interesting is the direct relation of the Pauli repulsion (ΔE Pauli) between the cations and the GQ scaffold, the strain energy (ΔE strain) required to accommodate these cations in their cavity, and the ionic radius of M (see Figure 7). We have deliberately excluded Li+, Mg2+, and Zn2+ from this analysis since the coplanar mode of their coordination leads to a non‐straightforward comparison of the values of ΔE Pauli. It is apparent how, with the growing of ionic radii and the steric repulsive forces inside the channel, the ΔE strain term decreases. In the previous section, we demonstrated that smaller cations (i. e., smaller ionic radius) lead to more shrinkage of the GQ cavity, and therefore to a higher strain energy (ΔE strain) than for the larger cations. The reason why smaller cations require more deformation lies in the fact that these lighter atoms have less diffuse orbitals to interact with the GQ scaffold and must approach the guanine bases much closer, which is reflected by the shorter d[O⋅⋅⋅M ] in Table 1. Less diffuse orbitals lead at the same time to a less destabilizing ΔE Pauli for the smaller metal cations, as the filled orbitals have less overlap with the electron density of the GQ cavity.

Figure 7

Relation between the cationic radii, strain energy of the GQ scaffold (blue lines), and Pauli Repulsion between the GQ cavity and cations (red lines). The dotted lines between the data points are there to guide the reader. Energies (in kcal mol−1) are calculated at ZORA‐BLYP‐D3(BJ)/TZ2P in the gas phase. The ionic radii (in Å) are taken from Ref. [21] and the data for GQ‐M1+ is taken from Ref. [8].

Kohn‐Sham Molecular Orbital Analysis of Model Systems

In this section, we examine the nature of the bonding between the metal cations and the GQ scaffold. Our Kohn‐Sham molecular orbital (MO) analysis provides insight into the origin of the more stabilizing orbital interactions (ΔE oi) for the alkaline earth metals and (post‐) transition metals within the quadruplex, compared to the alkali metals (Figure 6). Furthermore, we explain the observation that ΔE oi becomes less stabilizing when moving down Group 2 from Ca2+ to Sr2+, while an opposite trend is found for Group 12, where ΔE oi becomes more stabilizing from Cd2+ to Hg2+. In our MO analysis, we highlight how the s, p, and d lowest unoccupied molecular orbitals (LUMOs) on the metal cations are all involved in the bonding mechanism and contribute to the structure and stability of the GQ‐M complexes.

To understand the bonding mechanism between the cation and the GQ scaffold, a simplified model system was chosen by replacing the guanine residues in the scaffold of a GQ‐M (with M =Na+, K+, Rb+, Ca2+, Sr2+, Ba2+, Cd2+, Hg2+, and Pb2+) by formaldehyde molecules (see Figure 8 for a representative example). The O⋅⋅⋅M distances were fixed to be equal to the average distance in the original GQ‐M systems, while the formaldehyde molecules were constrained in the plane of the oxygen atoms (corresponding to S 8 symmetry) and then the system was optimized in C 2 symmetry at ZORA‐BLYP‐D3(BJ)/TZ2P. The symmetry reduction in the calculations is required as the higher symmetry S 8 is not implemented in the ADF program. The results of the energy decomposition analysis for the S 8 symmetric model system using C 2 symmetry are reported in SI Table S5 and are in line with the trends found for the GQ‐M systems in Figure 6.

Figure 8

Formaldehyde scaffold model system with S 8 symmetry (C 2 in computation) interacting with Ca2+.

To gain insight into the donor‐acceptor interactions between the cations and the GQ scaffold, the orbital energies of s, p, and d LUMOs of the cations along with the corresponding gross Mulliken populations were computed at ZORA‐BLYP‐D3(BJ)/DZ (Figure 9). The smaller DZ basis set was chosen over the larger TZ2P basis set as it is more suitable for the combination of compact lone pair orbitals of formaldehyde, and the more diffuse virtual orbitals of the metal cations, resulting in more reliable gross Mulliken populations that correlated with the EDA results (SI Table S5).

Figure 9

Metal cation M s, p, and d LUMO energy levels (in eV), and gross Mulliken populations (in electrons between brackets) resulting from the interaction with the central cavity in the S 8 symmetric formaldehyde scaffold model system, computed at ZORA‐BLYP‐D3(BJ)/DZ in the gas phase with C 2 symmetry. Note that with the DZ basis set, Na+ has no d unoccupied orbital.

From Figure 9 it becomes evident that the metal virtual levels are lowered in energy for the divalent alkaline earth metal cations and even further for the transition metal cations, compared to the monovalent alkali cations. This orbital lowering is in line with the increased ionic charge, and the increasing electronegativity along this trend. The low‐energy virtual orbitals of Group 2, 12, and 14 can accept more electron density from the GQ scaffold, reflected by the higher gross Mulliken populations, thus leading to more favorable orbital interactions compared to Group 1.

Figure 9 also shows that the s, p, and d unoccupied orbitals of Group 1 and 2 metal cations all engage in strong donor‐acceptor interactions with the occupied orbitals of the simplified formaldehyde scaffold, reflected by the gross Mulliken populations. This can be understood from the different in‐phase and out‐of‐phase combinations of the occupied orbitals on the formaldehyde scaffold, as is further highlighted in the next paragraph. Interestingly, only the s and p LUMOs of Cd2+ and Hg2+ of Group 12 are populated, as the unoccupied d orbitals are too high in energy and possess neglectable gross Mulliken populations. In the case of Pb2+ of Group 14, only the p level is receiving electron density from the formaldehyde scaffold, with the highly destabilized s and d orbitals playing no significant role in the coordination.

The highest occupied molecular orbitals (HOMOs) of a single layer of the formaldehyde model system are depicted in Figure 10. These HOMOs, with a large coefficient on the carbonylic oxygen atoms, make up the key electron‐donating orbitals to interact with the LUMOs on the metal cations due to their most pronounced overlap with the M virtual orbitals. The reason why the s, p, as well as the d unoccupied orbitals can engage in donor‐acceptor interactions (vide supra) is that the occupied orbitals of the formaldehyde monomers form different in‐phase and out‐of‐phase combinations, which allows the formation of both σ and π‐type bonds between the scaffold and the cation. For example, the HOMO−3 (Figure 10) has all the monomers’ filled orbitals in phase in the central cavity, which makes it a suitable bonding partner for the spherical s orbitals of the metals. This concept can be extended to the stacked quartets, which results in bonding and antibonding combinations that are either suitable for orbital interactions with s, p, or d metal LUMOs. So, the strength of the donor‐acceptor interactions depends on both the energy of the M LUMO orbitals, as well as their overlap with the quartet's filled orbitals.

Figure 10

HOMO orbitals (isosurfaces at 0.03 a.u.) of one layer of the S 8 symmetric formaldehyde scaffold model system, that are involved in metal cation coordination, computed at ZORA‐BLYP‐D3(BJ)/TZ2P in the gas phase with C 2 symmetry.

To understand which unoccupied orbitals of the metal cations are mainly responsible for the orbital interaction (ΔE oi) destabilization when moving from Ca2+ to Sr2+, while ΔE oi becomes more stabilizing going from Cd2+ to Hg2+, we examine the importance of the various metal virtual orbitals (s, p, and d) to the cation binding mechanism in terms of their contribution to ΔE oi. Again, we make use of the simplified formaldehyde model system, but in this case, we rotated one formaldehyde layer by 45° to obtain a C 4h instead of S 8 symmetric system (Figure 11). The sub‐symmetry C 2h was used in the calculations as the higher symmetry C 4h is not implemented in the ADF program. The advantage of using C 2h symmetry in the calculations is that the virtual s, p, and d orbitals of the metal belong to different irreducible representations, which allows us to separate the p orbital interactions from the s and d orbital interactions. In the formaldehyde model system in Figure 11, the O⋅⋅⋅M distances are fixed to the average distance in the original GQ‐M system. Furthermore, the formaldehyde monomers in the top and bottom layer are constrained to be in the same plane as the oxygen atoms and are directly on top of each other. This approach is justified, as the EDA results of the C 2h model systems produce similar energetic terms with identical trends as the C 2 model system. (SI Table S5).

Figure 11

Formaldehyde scaffold model system with C 4h symmetry (C 2h in computation) interacting with Ca2+.

The orbital interactions per irreducible representation (irrep) are gathered in Table 2. As the contribution from the A irrep is the only irreducible representation for Cd2+ to Hg2+ that follows the same trend as the total orbital interaction (ΔE oi), we can conclude that the trend in ΔE oi comes entirely from the A irreducible representation. So, the trend is determined by the orbital interactions involving either the metal s or d LUMOs, and not by the interactions with the p virtual orbitals. We cannot separate the s and d orbital interactions completely because the s, , , and orbitals all belong to the fully symmetrical A irrep. However, the empty d orbitals of the Group 12 transition metals Cd2+ and Hg2+ are too high in energy to accept significant electronic density (Figure 9) and are thus expected to not contribute to ΔE oi(A). Therefore, the stabilization of ΔE oi when going from Cd2+ to Hg2+ originates from differences in the virtual s orbitals. The ns LUMO is significantly stabilized from Cd2+ to Hg2+ (Figure 9), which can be ascribed to the significant increase of the nuclear charge by +32 a.u., which goes along with a large increase in the electrostatic as well as relativistic stabilization of the ns valence atomic orbital. On the other hand, the nuclear charge increases only by +18 a.u. from Ca2+ to Sr2+, leading to no net stabilization of the 5 s valence orbital of Sr2+ compared to the 4 s of Ca2+. The larger increase in nuclear charge from Cd2+ to Hg2+, compared to that from Ca2+ to Sr2+, results from the fact that Hg is the first Group 12 element that is preceded in its period by 14 f‐block elements (i. e., the lanthanides). Thus, the nature and energy of the metal's unoccupied atomic orbitals are of great relevance for determining the strength of the orbital interactions and therefore the interaction with the internal channel of guanine quadruplexes.

Table 2

Total orbital interaction (ΔE oi) and contributions (in kcal mol−1) from the orbitals of the A (s and d), B (d), A (p), and B (p) irreducible representations in the C 4h symmetric formaldehyde model systems interacting with Ca2+, Sr2+, Cd2+, and Hg2+ computed at ZORA‐BLYP‐D3(BJ)/TZ2P in the gas phase with C 2h symmetry.

	Ca2+	Sr2+	Cd2+	Hg2+
Ag (s, dxy , d_z^2 , d_x^2-y^2 )	−47.7	−42.4	−72.7	−102.8
Bg (dxz , dyz )	−42.5	−38.3	−23.4	−20.1
Au (pz )	−26.2	−24.1	−33.0	−30.3
Bu (px , py )	−31.2	−26.9	−45.6	−41.0
Total	−147.6	−131.7	−174.7	−194.2

A (s, d, , )

B (d, d)

A (p)

B (p, p)

Conclusions

Divalent pollution metal cations are proposed to induce genotoxicity by strong coordination to guanine quadruplex (GQ) DNA. The affinity of the internal channel site of GQs towards divalent cations (M2+) of the alkaline earth metals (Group 2: Mg2+, Ca2+, Sr2+, and Ba2+) and (post‐)transition metals (Group 12: Zn2+, Cd2+, Hg2+, and Group 14: Pb2+) was quantified using dispersion‐corrected density functional theory (DFT‐D) computations with implicit aqueous solvation. The results for the divalent cations are compared to the monovalent alkali metal cation series (Group 1: Li+, Na+, K+, and Rb+), which includes the cations Na+ and K+ that coordinate to quadruplexes under biological conditions.

The bond energy (i. e., cation affinity) results from a competition between a large stabilizing interaction energy and a large destabilizing desolvation energy. Although the strain energy varies relatively less among the different cations, smaller cations induce a higher degree of GQ cavity shrinkage upon coordination, leading to a higher strain energy. Interestingly, at the same time, compact elements experience less steric Pauli repulsion within the GQ cavity due to less diffuse filled orbitals. Furthermore, smaller cations encounter more destabilizing desolvation energies, but can in general engage in more favorable orbital interactions with the GQ cavity. Overall, this leads to an optimal bonding energy for K+ (Group 1), Sr2+ (Group 2), Hg2+ (Group 12), and Pb2+ (Group 14), of −53.2, −69.3, −71.1, and −74.2 kcal mol−1, respectively. The enhanced cation affinity for the divalent cations compared to K+ results primarily from more favorable electrostatic and orbital interactions. It follows from our Kohn‐Sham molecular orbital analysis that the energies of the unoccupied orbitals of the divalent cations M2+ are lower in energy compared to the monovalent alkali cations, which is in line with the increased ionic charge and electronegativity of these elements. These low‐energy unoccupied levels can participate in stronger donor‐acceptor interactions within the GQ cavity, leading to enhanced orbital interactions and thus higher affinities. The donor‐acceptor orbital interactions involve s, p, and d virtual orbitals for the Group 1 and 2 metal cations, s and p for Group 12, and only p orbitals for Group 14, while commonly only the metal s orbitals are considered relevant in the GQ‐M interaction.

Our results are in line with the experimental observation that divalent cations, and especially Hg2+ and Pb2+, induce GQ formation and coordinate with high affinities to the GQ cavity, possibly leading to the irreversible displacement of K+ and Na+ in vivo. This work contributes to a better understanding of how divalent metal cations could cause genomic instability after mercury and lead poisoning by coordination to GQ DNA.

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supporting Information

Click here for additional data file.