Low-frequency vibrational modes in G-quadruplexes reveal the mechanical properties of nucleic acids†

Low-frequency vibrations play an essential role in biomolecular processes involving DNA such as gene expression, charge transfer, drug intercalation, and DNA–protein recognition. However, understanding the vibrational basis of these mechanisms relies on theoretical models due to the lack of experimental evidence. Here we present the low-frequency vibrational spectra of G-quadruplexes (structures formed by four strands of DNA) and B-DNA characterized using femtosecond optical Kerr-effect spectroscopy. Contrary to expectation, we found that G-quadruplexes show several strongly underdamped delocalized phonon-like modes that have the potential to contribute to the biology of the DNA at the atomic level. In addition, G-quadruplexes present modes at a higher frequency than B-DNA demonstrating that changes in the stiffness of the molecule alter its gigahertz to terahertz vibrational profile.


SUPPLEMENTARY REFERENCES 43
Electronic Supplementary Material (ESI) for Physical Chemistry Chemical Physics. This journal is © the Owner Societies 2021 Figure S1. Optical Kerr-effect spectra of aqueous Na2(5'-GMP) solutions at different concentrations and temperatures showing changes in the spectra due to the formation of stacks and G-quadruplexes. Temperatures run from 283 (blue) to 328 K (red) in steps of 5 K and from 328 K to 358 K (orange) in steps of 10 K. The fitting deviations (       The temperature runs from 283 (blue) to 328 K (red) in steps of 5 K and from 328 K to 358 (orange) in steps of 10 K. Below the residuals between each spectrum and the fittings are shown. The remaining rows contain the spectrum for each case (grey) at 283 K (third row) and 358 K (fourth row) and the functions employed to fit the spectra. Diffusion associated functions (Debye and Cole-Cole) are coloured in red and vibration associated functions (Brownian oscillators) in orange. The blue band in the GMP·Li + and GMP·Cs + spectra at 283 K is the result of the sum of the two Brownian oscillators associated with GMP stack conformation (BA + BB). (Linear scale version of Figure 5).       Table S1. Fit parameters used to fit the high-frequency part of all the OKE spectra. Frequency (ω/2π) and damping rate (γ/2π, in brackets) measured in THz of the Brownian oscillators used to model the vibrational spectra of nucleotides and oligomers with G-quadruplex (G4), single-stranded (ss), and double-stranded (ds) conformations. 307.5 ± 0.8 Table S2. Melting points for the G4 conformation in Na2(5'-GMP) at different concentrations obtained for the curves plotted in Supplementary Figure 6.

Supplementary Note 1: Water OKE spectra and subtraction procedure
Measured OKE spectra have contributions from the solvent and the solvated nucleic acids. In order to extract the latter, the spectrum of water was directly subtracted from the measured OKE spectra (Supplementary Figure 11).
(2) To simplify the analysis, the spectra of water were simulated in the subtraction process using a model of the dependence of the water spectra on temperature. This model was created after the measurement and fitting of the OKE spectra of water between 283 and 368 K at intervals of 5 K (Supplementary Figure 12). Each spectrum was fitted using a combination of one Cole-Cole function for the diffusive translational motions of water molecules, two gaussian oscillators for the LA and TA phonon modes and one Brownian oscillator for the librational mode at 12.5 THz (Supplementary Table 1

⇌ $
We can estimate the enthalpy and entropy change during the process using van 't Hoff equation for the equilibrium constant: Substituting, we obtain: Formally, the denaturation temperature, or melting point (Tm) is defined as the temperature at which half of the G4 structure is present as free GMP. This is: Replacing at van 't Hoff equation:

Supplementary Note 3: Computational details
All computations were carried out in Gaussian 09 program. (6) PBE0 density functional (7, 8) with Grimme dispersion correction (9) combined with def2-SVP basis set (10) was used as a default method. PBE0 is a hybrid functional shown to provide a satisfactory description for spectroscopic properties (8), while def2-SVP basis set usually allows obtaining adequate estimations for geometry parameters for relatively large molecules consisting of the first and second-row elements in a substantially small CPU time. To check the structural accuracy of PBE0/def2-SVP, M06-2X (11) calculations in def2-TZVP (10) basis set were carried out. The solvent effects were accounted for by using the polarizable continuum model (12) (IEFPCM as implemented in Gaussian) with the dielectric constant εstatic = 78 (water). We do not aim to quantitatively reproduce OKE spectra in this work; instead, we aim to provide an insight into what type of nuclear motions could be responsible for high-intensity signals in the low-frequency part of Raman and OKE spectra. In order to do that, we constructed 1-layer, 2-layer and 3-layer models for G4 complexes with sodium and potassium (see Figure S18). For complexes with potassium, the cation was placed in-between the layers, while for complexes with sodium, the cation was in the plane of the G4 layer. Initial geometry optimizations were performed without symmetry; resulting structures were found to be close to S8 (2-layer systems), C4 (3-layer system with Na + ) and C4h (3-layer system with K + ) point groups and symmetrized using Chemcraft software (12) and optimized with symmetry constraints. Harmonic analysis on these structures showed that all Hessian eigenvalues are positive (i.e., there are no imaginary frequencies). The 3-layer model is arguably the smallest reasonable model to show most types of nuclear motion patterns expected in larger stacks (it has two "outer" and one "inner" layer). The majority of calculations was done in vacuum; in order to estimate the solvent effects on calculated spectra, PCM calculations for 2-layer systems were performed.
The OKE intensity for mode i was estimated as follows: is an anisotropic polarizability derivative along the mode i, and M ) stands for the mode wavenumber, ) is a depolarization ratio and ) is a Raman activity. Spectra were broadened with gaussians with 15 cm -1 half-width.

Structure
From the analysis of Table S3 and Table S4 one can see that the equilibrium geometry parameters of the guanine molecule obtained by PBE0-D3/def2-SVP and M06-2X/def2-TZVP are very similar. Moreover, changes in intramolecular distances and bond angles when comparing free guanine and G-quartet without cation are consistent for these two methods. Among the most notable structural changes for the guanine unit is the significant increase in the CO distance coupled with the decrease in adjacent CN distance, resembling a keto-enol tautomerism resulting from a rearrangement of the electron density in the vicinity of the oxygen atom ( Figure S15, the equilibrium is shifted to II in complexes). Thus, these geometry parameters could serve as local structural descriptors providing an additional argument to distinguish free guanine from that in complexes, e.g., in XRD data.
Non-bonding interactions in the G-quartet without the cation seem to be more difficult to reproduce: PBE0 with def2-SVP and def2-TZVP as well as M06-2X with def2-TZVP show higher discrepancies in hydrogen bond angles and distances. Moreover, while M06-2X provides a planar C4h structure for the G-quartet, a nonplanar S4 equilibrium geometry is obtained by PBE0-D3 both in def2-SVP and def2-TZVP basis sets. Looking at the H-N hydrogen bond length (Table S3 and Table S4), one could speculate that for def2-TZVP structures, the "outer" circle of hydrogen bonds (i.e., H16···N2') does not exist, unlike for def2-SVP ( Figure S16). This is indicated by notably larger distances (r((N6)H16···N2') is about 2.5 Å while the typical H···Y distance is 1.6-2.0 Å). Hence, among the three different methods used -PBE0-D3 with def2-SVP and def2-TZVP and M06-2X with def2-TZVP -each provides a different geometry, indicating that not only the DFT functional but also the basis set could lead to notably different structures. However, the energy difference between the C4h transition state in PBE0-D3/def2-TZVP (very similar to the minimum in M06-2X) and the nonplanar S4 minimum predicted by the same method is only 0.13 kcal/mol, which, given the accuracy of DFT methods, is essentially zero. Apparently, the potential energy surface section along the corresponding intermolecular mode is very shallow (i.e., large changes in geometry are accompanied by small changes in energy), and in that case one should look at the wave function (or probability density) of (at least) the zero vibrational level to resolve the planarity problem (12). Thus, to provide a better understanding of the intermolecular interactions in the Gquartet without a cation, a benchmark study would be desirable, which however goes beyond the scope of this paper.
The structure of the G-quartet with cation is much less sensitive and generally depends on the cation ( Figure S17). One can see that the G-quartet with K + has a conical shape while with Na + , the structure is very close to planar. These calculations already indicate where the cations are located in extended Gquadruplexes: K + is in-between layers, while Na + lies approximately in the layer (in agreement with experiment). (13) Since DFT functionals normally describe the interactions with a cation in a reasonable and robust way, we believe that PBE0-D3/def2-SVP can provide a qualitatively correct picture of the structure and vibrational properties of the G4 complexes with cations. Comparing hydrogen bond distances, one can see that they are quite different in G-quartets with and without cation, hinting at the governing role of cation in the layer geometry (i.e., relative positions of guanines within a layer).
Analysing changes in structural parameters for G-quadruplexes with alkali cations, it should be mentioned that the hydrogen bond distances for different layers are also different; quite expectedly, in 2layer complexes, they are closer to the "outer" layers in 3-layer complexes. One can also note, along with changes in hydrogen bonds geometry, changes in CO distance: the latter is larger in the outer layers of 3layer structures, and very similar to that in "free" G4 in the inner layer. Average interlayer distances without cations in Å (and distances between O4 planes, that could be thought of "minimal" interlayer distances) are as follows: 3.22 (2.95); 3.23 (3.07); 3.26 (3.04); 3.24, 3.25 (2.99, 3.06) for 2-and 3-layer complexes with K + and Na + , respectively (two values are given for 3-layer complex with Na + because of lower symmetry). These values should correspond to distances between cations in a sufficiently long stack; for our systems, K-K distance is 3.46 Å (from 3-layer complex), and Na-Na are equal to 3.34 (2-layer), and 3.30 and 3.02 Å (in 3-layer system). The calculated K-K distance is in good agreement with the reported experimental value for the 4-layer complex formed in DNA (3.5 Å) (14). For 3-layer complexes, it is to be expected that, due to guanine molecules being non-equivalent in outer and inner layers, notable splittings in vibrational spectra may occur.

Vibrational spectra
Indeed, one can notice that with adding a layer, low-frequency (LF) part of OKE spectra becomes more complicated ( Figure S19, Figure S21, Figure S23). A tentative assignment based on visual inspection of the vibrational motions reveals (see Table S9, Table S10, Table S11, and Table S12) that the Raman active modes can be divided into three groups: 1) < 40 cm -1 -interlayer modes (dispersion interaction between layers) 2) <150 cm -1 (excluding the previous group) -modes involving various displacements of guanine molecules as a whole near their respective equilibrium positions; these could also be described as collective hydrogen bond stretching / bending vibrations. 3) > 150 cm -1 -collective intramolecular modes could be approximately described as linear combinations of normal modes of a free guanine molecule. Modes belonging to the first two groups are purely intermolecular; as one can see, they clearly dominate in LF part for complexes with potassium, to a lesser extent -with sodium. For the first group, there is just one (degenerate) mode that could be easily described as an interlayer displacement. According to our calculations, this mode is present (in this context it means it shows a strong relative intensity) in all studied systems while having the highest intensity in all but 3-layer complex with sodium. Comparing 2-and 3-layer structures, it could be noted that corresponding harmonic frequency decreases (on 10 cm -1 ); also, it is notably lower in complexes with sodium. From a theoretical viewpoint, this kind of motion is important since it could lead to instability of a large enough stack of guanine layers (supposing there are no other constraints preserving the structure, e.g. side chains in nucleic acid). Indeed, calculations performed for a 2-layer stack with sodium in the PCM model with water (utilizing C4 symmetry constraint) resulted in a large imaginary frequency corresponding to this kind of motion. Conversely, the 2-layer complex with potassium optimized within the same procedure resulted in dynamically stable structure that is very similar to that in vacuum.
It is peculiar that intralayer "breathing" mode has shown a notable (still weak) signal only for the 2-layer complex with sodium at approximately 85 cm -1 .
There are quite a few modes falling into the second group with their number increasing while adding layers. They are for the most part similar, and the main difference arises from how the amplitudes of corresponding displacements are redistributed within and in between layers. Usually, they are delocalized for the whole structure; however, due to notable structural discrepancy between inner and outer layers, say, modes 17 and 19 in 3-layer complex with potassium are "layer-specific". One could speculate that the increase of the number of layers will lead to less structural deviation and hence, more delocalization for these types of nuclear motions. Having in mind the above discussion of the accuracy (regarding the description of hydrogen bonds) of the method we utilize we do not think that a more sophisticated assignment would be meaningful.
While it is not in the focus of this work, it might also be appropriate to briefly discuss types of nuclear motions leading to higher intensity peaks in the low-frequency part of IR spectra. Bands peaking at approximately 140 and 250 cm -1 for all studied G-quadruplexes could be assigned to (ungerade) collective G-G ipb and ν3 vibrations and are similar (albeit with a different symmetry) to those observed in OKE spectra. The lowest frequency vibration with notable IR intensity could be tentatively described as Tz mode for cations and their local surrounding and appears in 30-40 cm -1 region for complexes with K + and 50-60 cm -1 for complexes with Na + . Finally, a mode representing nuclear pattern that does not give rise to notable Raman activity, is a degenerate Tx, Ty mode for cations only. I.e., cation "string" is displaced relative to layers. This mode is interesting for not only being (almost uniquely for these systems) well localized but also because its position strongly depends on the cation: it is about 180 cm-1 for both complexes with K + , 258 cm -1 for 2-and 232 cm -1 for 3-layer complex with Na + . It is important to understand that relative intensities of low-frequency (< 350 cm -1 ) to high-frequency part for IR spectra are very small, so all discussed IR features may be very difficult to verify experimentally.
While it is possible to provide an intuitive explanation upon relative intensities in IR spectra, rationalizing Raman (or OKE) spectra is a less straightforward problem. In order to provide at least basic considerations, one can look at the symmetry of corresponding modes. In dipole approximation, fundamental transitions belonging to the following representations are allowed in Raman: 1) S8 point group (2-layer systems): A, E2, E3.
3) C4 (3-layer with Na + ): A, B, E (all). Analysing Table S9, Table S10, Table S11, and Table S12, it is evident that for all complexes normal modes belonging to degenerate representations (E3, Eg, E) lead to higher relative intensities in the Raman spectrum. However, many strong signals in the computed spectrum of the 2-layer complex with Na + belong to a fully symmetric (A) representation. Thus, it is easy to show that most Raman-active modes have a rotation symmetry: (Rx, Ry) transform as E3, Eg, E in S8, C4h, and C4 point groups, respectively, while Rz transforms according to A (or Ag) irreducible representation. The fact why the 2-layer complex with sodium shows somewhat different behaviour remains unexplained.
Finally, comparing calculated OKE spectra for G4 complexes in vacuum and in polarized media (Figure S19 and Figure S23) one can see that solvent effects should definitely be accounted for when aiming to reproduce the intensity pattern in LF part, e.g., the relative intensity of interlayer displacement mode is smaller in water than in vacuum. However, structure, vibrational patterns and harmonic frequencies are very similar for these two complexes; thus, we believe that our qualitative conclusions would in general be supported by future in-depth studies.
In order to provide a more rigorous, quantitative theoretical description for the systems in question, one should 1) use a better method (def2-SVP basis set is too small to provide a good description of a hydrogen bond and polarizability); 2) investigate larger stacks (probably using PBC calculations, or better, modelling finite stacks with a number of layers close to observed in our experiment); 3) go beyond double harmonic approximation when it comes to spectra. Figure S14. Atomic enumeration in guanine fragment.