The benchmark of 31P NMR parameters in phosphate: a case study on structurally constrained and flexible phosphate

A benchmark for structural interpretation of the P NMR shift and the JP,C NMR spin–spin coupling in the phosphate group was obtained by means of theoretical calculations and NMR measurements in diethylphosphate (DEP) and 5,5-dimethyl-2-hydroxy-1,3,2-dioxaphosphinane 2-oxide (cDEP). The NMR parameters were calculated employing the B3LYP, BP86, BPW91, M06-2X, PBE0, KT2, KT3, MP2, and HF methods, and the 6-31+G(d), Iglo-n (n = II, III), cc-pVnZ (n = D, T, Q, 5), aug-cc-pVnZ (n = D, T and Q), and pcS-n and pcJ-n (n = 1, 2, 3, 4) bases, including the solvent effects described with explicit water molecules and/or the implicit Polarizable Continuum Model (PCM). The effect of molecular dynamics (MD) on NMR parameters was MD-calculated using the GAFF force field inclusive of explicit hydration with TIP3P water molecules. Both the optimal geometries and the dynamic behaviors of the DEP and cDEP phosphates differed notably, which allowed a reliable theoretical benchmark of the P NMR parameters for highly flexible and structurally constrained phosphate in a one-to-one relationship with the corresponding experiment. The calculated P NMR shifts were referenced employing three different NMR reference schemes to highlight the effect of the P NMR reference on the accuracy of the calculated P NMR shift. The relative Dd(P) NMR shift calculated employing the MD/B3LYP/Iglo-III/PCM method differed from the experiment by 0.16 ppm while the NMR shifts referenced to H3PO4 and/or PH3 deviated from the experiment notably more, which illustrated the superior applicability of the relative NMR reference scheme. The JP,C coupling in DEP and cDEP calculated employing the MD/B3LYP/ Iglo-III(DSO,PSO,SD)/cc-PV5Z(FC)/PCM method inclusive of correction due to explicit hydration differed from the experiment by 0.32 Hz and 0.15 Hz, respectively. The NMR calculations demonstrated that reliable structural interpretation of the P NMR parameters in phosphate must involve both the structural and the dynamical components.


Introduction
The phosphate group is the key element within the sugarphosphate backbone of nucleic acids. Although the major stabilizing effects in nucleic acids come from hydrogen bonding and stacking interactions between bases and base-pairs, their 3D framework is given by the conformation of the backbone. [1][2][3][4] Orientations of the torsion angles adjacent to the backbone phosphate are distinctive for backbone conformational classes in nucleic acids. 5 Much of the backbone flexibility originates from rotations around the phosphodiester linkage and that is why interconversions among distinct conformations of the backbone progress typically via the phosphodiester group. Structural resolution of the backbone phosphate is therefore central to all structural studies concerning nucleic acids.
To date, the majority of the known structures of nucleic acids have been resolved in crystals by employing X-ray crystallography methods. In liquids, NMR spectroscopy is employed instead. The NMR parameters including the Nuclear Overhauser Effect (NOE), NMR chemical shift of atom A (d(A)), indirect NMR spin-spin coupling between atoms A and B separated by n bonds ( n J A,B ), and residual dipolar coupling (RDC) can provide structural constraints for obtaining the 3D structure of DNA and RNA molecules. 6 Moreover, probably the only method for the reliable refinement and testing of MD force fields can be currently done with structural data derived from NMR experiments. 7 Hence, knowledge of the rules for the structural interpretation of NMR parameters is indispensable. The empirically derived rules can be employed; however, their accuracy need not be guaranteed in contrast to the ab initio calculations of NMR shifts and J-couplings. 8 As for the nucleic acids, the NMR measurements and calculations highlighted the base pairing, the orientation of the base with respect to the sugar, the structure of the sugar-phosphate/phosphonate backbone in canonical/ modified oligonucleotides aimed at medicinal applications, the oxidative damage of DNA, and the 3D structure of metallo-DNAs including Hg/Ag-mediated base pairs. [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25] The 31 P NMR measurements conveniently employ the 100% natural abundance of the isotope. Recent calculations of the 31 P NMR shifts in organophosphorus compounds differing by hundreds of ppm illustrated the reliability of the theoretical DFT methods. 26 By contrast, the 31 P NMR shifts in the nucleic acids differ typically at most by a few ppm which makes their assignment and interpretation truly challenging. [27][28][29][30] The impact of the accurately interpreted 31 P NMR parameters on the quality of NMR-resolved structures of nucleic acids is significant. The theoretical 31 P NMR parameters were therefore intensively studied. Already pioneering calculations have unveiled the dependence of the 31 P NMR shift on the geometry of phosphate. 31,32 Later, by using advanced DFT calculations, Munzarová et al. illustrated that the A-RNA, BI/BII DNA, and Z-DNA backbone conformations are distinguishable with 31 P shifts. 33,34 Recently, Munzarová et al. calculated the dynamically averaged 31 P NMR shifts in RNA and B-DNA by employing the combined MD/DFT calculation method. 35,36 The calculations of 31 P NMR shift and 2 J P,C NMR spin-spin coupling in phosphate illustrated that various conformations of the backbone in nucleic acids can be distinguished with these NMR parameters. 37 The effect of solvation due to Mg 2+ on the NMR parameters was calculated. 38 These NMR calculations demonstrated that the inclusion of the solvent and the dynamics involved affect the magnitude of the NMR parameters in phosphate. In the relation of the theoretical 31 P NMR shift to experiments another issue should be considered; the effect of the 31 P NMR reference. The d( 31 P) NMR shift in nucleic acids is usually measured with respect to phosphoric acid; however, calculation of the d( 31 P) may lead to erroneous assignment and misinterpretation. Another 31 P NMR reference scheme, the secondary PH 3 reference, has nowadays become popular in theoretical 31 P NMR calculations. 39 General aspects of the 31 P NMR referencing were addressed in this work to unambiguously illustrate the performances of the currently available reference schemes and to demonstrate the effect on the structural interpretation of the 31 P NMR shift.
As for J-coupling, the NMR parameter provides the so called long-distance structural restraint due to its dependence on the geometries of all the atoms within the NMR spin-spin pathway. 40 In nucleic acids, the J-couplings were employed for NMR measurement/detection of the H-bond and for determination of the backbone conformations. [41][42][43][44][45][46] A closer look at the backbone phosphate unveils the lack of suitable atoms measurable by means of NMR, except for the phosphorus and the nearest carbon and hydrogen atoms. The 31 P NMR shift and 2 J P,C NMR coupling are therefore unique NMR parameters in the phosphate.
Theoretical methods for calculations of the 31 P NMR parameters have not been yet benchmarked because of the following reasons and obstacles. The previous 31 P NMR calculations aimed at nucleic acids employed simplified models of the phosphate group to include all necessary effects at the level corresponding at least to DFT. The plausibility of the di-methyl phosphate or methyl-ethyl-phosphate that were previously employed to mimic the phosphate group in nucleic acids is therefore unclear. On the other hand, high-accuracy NMR calculations in the nucleic acids are hardly doable even nowadays and the usage of high-precision methods is usually balanced with the adequate size of the molecular model. The DEP and cDEP molecules depicted in Fig. 1 include presumably the smallest possible model of the phosphate group, which is still chemically similar to the phosphate in nucleic acids. As can be anticipated from their chemical structures, the geometries and dynamics of the DEP and cDEP phosphates should allow controllable examination of the effects on 31 P NMR parameters. The benchmarked 31 P NMR parameters in DEP and cDEP therefore illustrate the accuracy of the calculation methods that are applicable for 31   2.1.3. The measurements of NMR spectra. The NMR spectra of DEP and cDEP were measured on a Bruker AVANCE-600 ( 1 H at 600.13 MHz and 13 C at 150.9 MHz) using a 5 mm TXI cryoprobe and a Bruker AVANCE-500 instrument ( 31 P at 202.3 MHz) using a 5 mm TBO BB-probe. About 10 mg of the sample was dissolved in 0.6 ml D 2 O and the NMR spectra were measured first at natural acidic pH and then as a strongly alkaline solution obtained by the addition of NaOD. The titration experiments showed the pK values of the phosphate groups in DEP and cDEP to be about 2. A trace amount of dioxane was added as the internal reference for the 1 H and 13 C spectra and the chemical shifts were recalculated to the TMS using d( 1 H) = 3.75 and d( 13 C) = 69.3 ppm. The 31 P NMR spectra were referenced to the external H 3 PO 4 (in capillary). The coupling constants J H,H and J H,P were determined from the 1D proton spectra. The carbon-13 chemical shifts and J C,P values were obtained from proton broadband decoupled 13 C spectra. The proton broad-band decoupled 31 P NMR spectra were used to get 31 P chemical shifts. The multiplets observed in gated proton-decoupled 31 P NMR spectra confirmed the J P,H values. The accuracies of the experimental chemical shifts and coupling constants are 0.01 ppm and 0.1 Hz, respectively.

The theoretical calculations
The geometry optimizations were carried out with the B3LYP 47-50 method, the 6-31+G(d) [51][52][53][54][55] atomic basis, and the polarizable continuum implicit solvation model (PCM) 56 describing the water solvent. The dependence of the energy on the orientation of the torsion angles adjacent to the phosphorus atom in the phosphate group was calculated with the B3LYP method and the 6-31G(d,p) basis by employing a 301 Â 301 geometric grid. The validity of all energy minima was checked with the vibration frequency calculation.
The molecular dynamic (MD) calculations employed the GAFF force field 57,58 and the ESP atomic charges 59,60 calculated using the HF method and the 6-31G(d,p) basis. The 50 Å Â 50 Å Â 50 Å cubic box in MD calculations included the DEP or cDEP molecule, the TIP3P 61 water molecules and one Na + ion to maintain the electroneutral state of the system. The system was first equilibrated within the two-step procedure including heating and 30 ps of free NPT MD calculation. The production NVT MD calculation included 35 ns under the standard laboratory conditions; the temperature was 300 K and the pressure was 1 atm. The MD snapshot geometries equally separated by 100 ps including water molecules within the first solvation shell were employed. The 350 MD snapshots were geometry optimized keeping the T1, T2, T3, and T4 torsion angles and geometries of the oxygen atoms of the water molecules fixed (Fig. 1). The optimized MD snapshots were used for dynamical averaging of the NMR parameters. The MD snapshots of the NMR reference molecules H 3 PO 4 and PH 3 were calculated similarly. The 400 snapshots for H 3 PO 4 and PH 3 neglecting explicit water molecules were geometry optimized keeping the O-P-O and H-P-H valence angles fixed. The enhanced MD sampling for the cDEP within 20 ns MD calculation at the temperature 400 K was used to study the interconversions between the conformers of energy minima.
The calculated 31 P NMR parameters included the 31 P NMR shift of the phosphorus atom and the 2 J P,C two-bond J-coupling between the phosphorus and carbon atom. The 31 P NMR shift was calculated employing the GIAO approach within the Coupled Perturbed (CP) DFT method. [62][63][64] The 2 J P,C was calculated employing the CP DFT method as a sum of the Diamagnetic-spin orbit (DSO), Paramagnetic-spin orbit (PSO), Fermi-contact (FC) and Spin-dipolar (SD) contributions. 65,66 The NMR parameters were calculated with the B3LYP method, the Iglo-III basis, 67 and the PCM water solvent. Only the FC term of 2 J P,C was calculated with the cc-pV5Z basis. 68 77,78 and HF 79 methods. The explicit water molecules in NMR calculations were described with the 6-31+G(d) basis. The effect of atomic basis on NMR parameters was calculated with the 6-31+G(d), Iglo-II, Iglo-III and cc-pVnZ (n = D, T, Q and 5), aug-cc-pVnZ (n = D, T and Q), and Jensen's pcS-n and pcJ-n, n = 1, 2, 3, 4 bases. [51][52][53][54][55][67][68][69][80][81][82][83] The NMR parameters calculated for the geometry of the global energy minima are called the static NMR parameters in the text. Accordingly, the statistically averaged NMR parameters are called the MD-averaged NMR parameters. The standard mean deviation from the statistical average S M was calculated for the MD-averaged NMR parameters as a sum where A i is the value of the respective NMR parameter in the i-th snapshot, and A av = P (A i /N) is the statistical average for N MD snapshots. The 31 P NMR shift was calculated using three different NMR references.
The d( 31 P) referenced to H 3 PO 4 was calculated as d( 31 P) = s( 31 P) (in H 3 PO 4 ) À s( 31 P) (in DEP or cDEP) where the s( 31 P) was the calculated 31 P NMR shielding. The d( 31 P) PH 3 using the secondary PH 3 reference was calculated according to van Wüllen as d( 31 P) PH 3 = s( 31 P) (in PH 3 , neglecting solvent) À s( 31 P) (in DEP or cDEP) À 266.1. 39 The relative Dd( 31 P) NMR shift was calculated as Dd( 31 P) = s( 31 P) (in DEP) À s( 31 P) (in cDEP) = d( 31 P) (in cDEP) À d( 31 P) (in DEP), which allowed unbiased comparison of the calculated 31 P NMR shift with the experiment neglecting the external NMR reference. The geometry optimizations and NMR calculations were carried out with Gaussian 09.D.01. 84 Only the KT2 and KT3 31 P NMR shifts were calculated with the Dalton 2016.0 program. 85 The MD calculations were carried out with Amber 10. 86

The geometry of the DEP and cDEP phosphates
The DEP and cDEP phosphates are chemically similar due to the equivalent chemical bonding in the vicinity of the phosphorus atom. The magnitudes of the 31 P NMR parameters therefore depend on the geometries of the DEP and cDEP phosphates. Two energy-equivalent minima were calculated for both DEP and cDEP because of the molecular symmetry (Table 1). In DEP, the orientations of the T2 and T3 torsion angles were both 72.81 and for the symmetry-associated energy minimum it was À72.81. We used the À1801 to 1801 convention rather than 01 to 3601 to display clearly the dependences of the energy and NMR parameters on T2 and T3. In cDEP, T2 and T3 were À59.21 and 59.21, and for the symmetry-associated minimum they were 59.21 and À59.21, respectively. The optimal geometries of DEP and cDEP thus differed notably.
The dependences of energy on T2 and T3 calculated for DEP and cDEP unveiled the geometries of the local energy minima (Fig. 2).
For DEP, the four symmetrically placed energy-minima were 0.8 kcal mol À1 above the global energy minimum (Fig. 2b). For cDEP, the two energy-minima were 3.8 kcal mol À1 above the global energy minimum (Fig. 2d). Importantly, the potential energy surfaces unveiled the rotational freedom of T2 and T3 in DEP and cDEP. The notable flexibility of the DEP phosphate that can be assumed from the sizable low-energy area within the potential energy surface was actually also MD-calculated (Fig. 2a). a The torsion angles in degrees calculated with the B3LYP method and 6-31G(d,p) basis (Fig. 1). b The energy minima for DEP and cDEP are depicted in Fig. 2b and d, respectively. The parameters for energy-equivalent symmetrical minima are in parentheses. The MD snapshots even reflected the curvature of the lowenergy profile, which illustrated the reliable MD sampling. By contrast, the flexibility of the cDEP phosphate was notably confined. The T2 and T3 rotational motions were constrained within a small area including the global energy minimum; elsewhere a sharp increase in energy was calculated (Fig. 2c). Interconversions between the two energy-equivalent energy minima of cDEP separated by the energy barrier of ca. 9 kcal mol À1 were not MD-calculated (Fig. 2c). The assumption of blocked interconversions, however, contrasted with the NMR experiment. The 13 C NMR shift of the methyl group and the 1 H NMR shift of the -CH 2 -group measured in cDEP clearly demonstrated fast interconversions between the energy minima (Table S1, ESI †). The energy barrier depicted in Fig. 2d thus can be affected due to the neglected involvement of other inner-ring torsional motions. The enhanced MD-sampling of the cDEP geometries calculated at a temperature of 400 K unveiled one accomplished interconversion ( Fig. S1 and S2, ESI †). The relative energies of the respective MD snapshots up to 10.7 kcal mol À1 illustrated the feasibility of the interconversions (Fig. S3, ESI †). The flexibility of the cDEP phosphate was notably constrained due to the molecular junction within the ring; however, fast interconversions between the two NMR-equivalent conformers occur.
DEP and cDEP phosphates differed notably as regards their optimal geometries and dynamic behaviour. None of the six DEP energy-minima geometries were similar to any energyminima geometries of cDEP. The flexibility of the DEP phosphate contrasted with the rigidity of the cDEP phosphate. The effect of the PCM implicit water solvent on the optimized geometries of DEP and cDEP was small. The calculated T2 and T3 torsion angles changed upon PCM hydration by less than 6.01 (Tables 1 and 2). The energy-optimized and MDcalculated geometries of cDEP were practically the same. For DEP, the MD-calculated average geometry was meaningless as the T2 and T3 averaged to E1801, whereas these orientations of T2 and T3 were in fact energy-forbidden. The cDEP and DEP phosphates thus can be regarded as geometrical antipodes; the calculated geometries of cDEP were not allowed for DEP and vice versa. The MD snapshot geometries describing the dynamic behaviour of DEP and cDEP were used for MD-averaging of the 31 P NMR parameters.

Preliminary NMR calculations -the basis set effect
The basis set dependence of NMR parameters was calculated using the 6-31+G(d), Iglo-n (n = II, III), cc-pVnZ (n = D, T, Q and 5), aug-cc-pVnZ (n = D, T and Q), and pcS-n and pcJ-n (n = 1, 2, 3, 4) atomic bases employing the B3LYP method and the PCM water solvent. The geometries of DEP and cDEP were optimized using the B3LYP method, 6-31+G(d) basis and the PCM water solvent if not written otherwise.
First, the s( 31 P) NMR shielding and 2 J P,C coupling in DEP and cDEP were calculated with Iglo-III for the grid point geometries used in the calculations of the potential energy surfaces in Fig. 2  (Fig. S13 and S14, ESI †). The s( 31 P) calculated for the global energy minima of DEP and cDEP ranged from 279.48 ppm to 383.37 ppm (Tables S2-S5, ESI †). The relative Dd( 31 P) NMR shifts ranged from À9.16 ppm to À0.67 ppm, except for the aug-cc-pVQZ where Dd( 31 P) was 5.52 ppm (Table S6, ESI †). The 2 J P,C coupling in DEP and cDEP ranged from À8.04 Hz to À5.13 Hz and from À7.13 Hz to À4.47 Hz, respectively.
The convergent-like behaviour of the NMR parameters was calculated for almost all basis sets; however, some NMR parameters obtained with the high-quality bases within particular series differed notably (Fig. 3). The relative Dd( 31 P) NMR shift within the pcS-n series converged rapidly, which illustrated the good performance of the Jensen's basis set. The basis set limit for Dd( 31 P) in this work with the pcS-4 was À2.14 ppm. The Dd( 31 P) within the cc-pVnZ series varied more and the Dd( 31 P) with the cc-PV5Z was À2.18 ppm. The Dd( 31 P) with Iglo-III was À2.17 ppm. The aug-cc-pVnZ set was not applicable for Dd( 31 P) calculation (Fig. 3a). The 2 J P,C coupling with the pcJ-n bases was calculated as the FC contribution with the pcJ-n plus the DSO, PSO and SD contributions with Iglo-III. Calculations of the three contributions with high-order pcJ-n bases would be impractical. The 2 J P,C in DEP and cDEP with the pcJ-n bases converged to À7.30 Hz and À6.56 Hz, respectively. The 2 J P,C in DEP and cDEP with the cc-pVnZ bases decreased to À7.02 Hz and À6.32 Hz, respectively. The Iglo-calculated 2 J P,C in DEP and cDEP increased to À7.63 Hz and À6.83 Hz, respectively. The 2 J P,C with aug-cc-pVQZ in DEP and cDEP was À6.62 Hz and À5.96 Hz, respectively.
The effect of geometry optimization on the NMR parameters was calculated by employing the atomic bases used in NMR calculations. The s( 31 P) NMR shielding in DEP and cDEP increased upon geometry reoptimization by up to 11.8 ppm, but for cc-pVDZ and aug-cc-pVDZ s( 31 P) decreased (Tables S3  and S5, Fig. S4, ESI †). The effect on the relative Dd( 31 P) shift was far smaller as the s( 31 P) in DEP and cDEP increased or decreased simultaneously (Fig. 3a). The Dd( 31 P) changed due to geometry reoptimization by less than 0.7 ppm and for cc-pV5Z it changed only by 0.1 ppm (Table S6, ESI †). The 2 J P,C mostly decreased, except for the aug-cc-pVDZ (Tables S3 and S5, ESI †). The effect of geometry optimization on 2 J P,C in DEP and cDEP ranged from À0.2 Hz to À0.7 Hz and from À0.1 Hz to À0.5 Hz, respectively. The basis set effect on MD-averaged NMR parameters was estimated for the more flexible DEP phosphate using only the first 20 MD snapshots (Fig. S5 and Table S7, ESI †). The s( 31 P) decreased upon MD-averaging, by up to 2.5 ppm with cc-pVTZ. The 2 J P,C increased, by up to 1.4 Hz with Iglo-II. The 2 J P,C calculated with cc-pV5Z increased upon MD-averaging from À7.0 Hz to À6.1 Hz. The effect of molecular dynamics on the NMR parameters should be therefore included.
The preliminary NMR calculations indicated suitable atomic bases. As regards Dd( 31 P), the high-quality bases within the conceivable series, pcS-4, cc-pV5Z and Iglo-III gave the same magnitudes of the relative NMR shift (Table S6, ESI †).
Out of the three bases, Iglo-III is well-suited for MD averaging of the 31 P NMR shift due to high accuracy for low computational demands. Regarding 2 J P,C , the basis set effect was examined for all of the four J-coupling contributions.
The FC contribution was prevalent and depended on atomic bases much more than the other three contributions where the DSO and SD were negligible and the PSO was de facto basis set independent (Tables S2 and S4, ESI †). The FC was determinative for the calculated magnitude of the 2 J P,C . To maintain high accuracy of MD-averaged 2 J P,C , the FC term was calculated with the cc-pV5Z basis while the remaining terms were calculated with the Iglo-III basis.

The 31 P NMR shift
The d( 31 P) NMR shifts measured in DEP and cDEP were referenced to H 3 PO 4 . In the theoretical calculations, three NMR references were employed. The d( 31 P) referenced to H 3 PO 4 was equivalent to the reference in the experiment. The d( 31 P) PH 3 was referenced to secondary NMR reference PH 3 according to van Wüllen. 39 Lastly, the relative Dd( 31 P) NMR shift was calculated as the mutual NMR shift of the two phosphates neglecting the external NMR reference. In this way, the effect of the NMR reference on the accuracy of the calculated 31 P NMR shift was unveiled.
The d( 31 P) NMR shift is usually measured in nucleic acids. The accuracy of the calculated d( 31 P), however, may be affected when calculation errors for s( 31 P) in the phosphate and in H 3 PO 4 are different. The effect can be anticipated due to different chemical bonding of the phosphorus atoms in phosphate and in H 3 PO 4 . By contrast, systematic errors in the calculations of s( 31 P) for chemically equivalent phosphates should likely cancel out in the calculation of the Dd( 31 P) NMR shift. The secondary PH 3 NMR reference was introduced to bypass the problematic H 3 PO 4 reference in theoretical calculations. 39 However, the calculation of small up-field and down-field 31 P NMR shifts such as in the phosphate still remained challenging. In the following, we will show that the choice of NMR reference may become critical for the assignment and interpretation of the 31 P NMR shifts.
The d( 31 P) measured in DEP and cDEP was 1.34 ppm and À2.68 ppm, respectively. The d( 31 P) calculated in DEP and cDEP with the B3LYP method and the pcS-4 basis was 2.14 ppm and 0.00 ppm, respectively. The Iglo-III-calculated d( 31 P) in DEP and cDEP was 2.03 ppm and À0.14 ppm, respectively. The cc-PV5Zcalculated d( 31 P) in DEP and cDEP was 2.80 ppm and 0.62 ppm, respectively. The calculated d( 31 P) shifts depended on highquality atomic bases rather notably. The d( 31 P) PH 3 shifts varied similarly; however, their magnitudes deviated from the experiment far more than those of the d( 31 P) shifts (Table 3). By contrast, the Dd( 31 P) calculated with the three high-quality atomic bases differed by less than 0.1 ppm, which indicated cancellation of the systematic errors in the calculated s( 31 P) for DEP and cDEP phosphates. The concurrent similar variations of the s( 31 P) in DEP and cDEP due to different atomic bases almost cancelled out, which made the Dd( 31 P) numerically stable (Table S6, ESI †). And vice versa, the s( 31 P) in H 3 PO 4 and PH 3 depended on atomic bases less systematically which explained the notable variations of the d( 31 P) and d( 31 P) PH 3 shifts (Table S8, ESI †). Importantly, the Dd( 31 P) calculated with the pcS-4, cc-pV5Z and Iglo-III bases demonstrated that accurate dynamical averaging of the Dd( 31 P) shift is achievable at low cost with the Iglo-III. The Dd( 31 P) calculated within the static approach involving the global energy minima of DEP and cDEP was À2.18 ppm whereas the measured Dd( 31 P) was À4.02 ppm.
The effect of MD averaging on s( 31 P) NMR shielding was noticeable (Fig. S6-S8 (Table S8, ESI †). The extent of fluctuations was reflected Fig. 3 The dependence of relative Dd( 31 P) NMR shift and 2 J P,C spin-spin coupling on the number of basis functions calculated with the 6-31+G(d), Iglo-II, Iglo-III, cc-pVnZ (n = D, T, Q and 5), aug-cc-pVnZ (n = D, T and Q), and pcS-n and pcJ-n (n = 1, 2, 3, 4) atomic bases for the geometries of DEP and cDEP molecules optimized with the 6-31+G(d) basis (filled symbols) and other atomic bases (open symbols). For the pcJ-n series, only the FC term of 2 J P,C was calculated with the bases, and the DSO, PSO, and SD terms were calculated with Iglo III. by statistical S M deviations. The S M for s( 31 P) in PH 3 was 0.9 ppm, while the S M for the other molecules did not exceed 0.2 ppm. Therefore, the d( 31 P) PH 3 changed upon MD-averaging more than the d( 31 P). The Iglo-III-calculated d( 31 P) in DEP and cDEP decreased upon MD-averaging by 1.17 ppm and 2.86 ppm, respectively. The MD-averaged d( 31 P) in DEP and cDEP was smaller than those in the experiment by 0.48 ppm and 0.32 ppm, respectively. When considering the error-bars corresponding to the sum of S M deviations for the respective s( 31 P), 0.23 ppm, the MD-averaged d( 31 P) agreed with the experiment. The MD-averaged d( 31 P) PH 3 in DEP and cDEP was larger than that in the experiment by 0.62 ppm and 0.78 ppm, respectively. The MD-averaged Dd( 31 P) was larger than that in the experiment only by 0.1 ppm. Such accuracy is acceptable for the reliable assignment and interpretation of the 31 P NMR shifts in nucleic acids. The effect of explicit water on MD-averaged Dd( 31 P) calculated with the 6-31+G(d) basis was only 0.1 ppm ( Table 3). The high accuracy and small dependence on both the atomic basis and the model of solvent make the relative NMR referencing superior as compared to the other two NMR reference schemes.
The geometry optimization of MD snapshots was necessary for achieving a high accuracy of MD-averaged 31 P NMR shifts. The usage of the original MD snapshots affected the magnitudes of the calculated s( 31 P) notably and their large variations were reflected by S M deviations (Table S8, ESI †). Consequently, the d( 31 P) and d( 31 P) PH 3 averaged with non-relaxed MD snapshots deviated from the experiment far more than their static and/or MD-averaged counterparts (Table 3). Strikingly, the Dd( 31 P) averaged with non-relaxed snapshots was only 0.9 ppm larger than that of the experiment and the deviation corresponded to the S M deviation. This illustrated again the favourable effect of the relative NMR referencing where the systematic geometrical effect on the NMR parameters due to the original geometry of the MD snapshots almost cancelled out.
The d( 31 P) in cDEP decreased upon MD-averaging by 2.86 ppm, whereas the d( 31 P) in the more flexible DEP decreased only by 1.17 ppm. This seeming contradiction can be explained by the dynamic effect on s( 31 P) (Table S8, ESI †).
The s( 31 P) in DEP and cDEP decreased upon MD-averaging by 2.13 ppm and 0.44 ppm, respectively. The larger dynamic effect on the local NMR property was thus actually calculated for the more flexible phosphate. The dynamic effect on d( 31 P) was due to the effect on the H 3 PO 4 reference where the s( 31 P) decreased upon MD-averaging by 3.3 ppm. The Dd( 31 P) can be interpreted more straightforwardly in terms of the geometry and dynamics of phosphate. The change of phosphate geometry when going from cDEP to DEP accounted for the ''static'' Dd( 31 P) À2.17 ppm ( Table 2). The increased flexibility of the phosphate when going from cDEP to DEP accounted for the ''dynamical'' Dd( 31 P) À1.69 ppm (Fig. 2). The dynamic effect on Dd( 31 P) can be highlighted in detail by considering the geometrical dependence of s( 31 P) (Fig. S13, ESI †). The s( 31 P) in the DEP MD snapshots was mostly below s( 31 P) for the energy minimum and s( 31 P) decreased upon MD-averaging (Table S8, ESI †). The MD geometries of cDEP were confined within a very narrow elliptic-like area where the main axis was perpendicular with respect to the s( 31 P) contour lines ( Fig. 2 and Fig. S13, ESI †), which resulted in a notable compensation of the dynamical variations of s( 31 P). Consequently, the s( 31 P) in cDEP decreased upon MD-averaging only by 0.44 ppm. The structure and dynamics of phosphate were elucidated with the theoretical calculation of the Dd( 31 P) NMR shift (Fig. 4).
The plausibility of the B3LYP-calculated 31 P NMR shifts was further supported by NMR calculations with other methods inclusive of PCM water hydration and the Iglo-III basis. The geometry optimizations were carried out with the method used in NMR calculation, the 6-31+G(d) basis and the PCM water implicit solvent. The MD-averaged Dd( 31 P) calculated with the BP86 method was À4.83 ppm and the static Dd( 31 P) was À2.88 ppm. The MD-averaged d( 31 P) in DEP and cDEP with BP86 was À1.15 ppm and À5.98 ppm, respectively. The B3LYP method thus provided a better agreement of the calculated Dd( 31 P) with the experiment than the BP86 method. The BP86 calculations, however, confirmed the superior performance of the relative NMR reference scheme. Moreover, the dynamic effects on Dd( 31 P) À1.95 ppm and À1.69 ppm were calculated  87 The KT2 calculations with the pcS-n bases converged rapidly (Table S10, ESI †). For increasing order of pcS-n, Dd( 31 P) was À3.11 ppm, À2.82 ppm, À2.99 ppm and À2.98 ppm. The KT2-calculated Dd( 31 P) with Iglo-III was À2.98 ppm, which demonstrated that the performance of Iglo-III is similar to the performance of pcS-4. The Dd( 31 P) values with the KT3 method were a little larger (Table S11, ESI †). For increasing order of pcS-n, the Dd( 31 P) was À2.65 ppm, À2.40 ppm and À2.59 ppm. The KT3-calculated Dd( 31 P) with Iglo-III was À2.58 ppm. Assuming the estimate of the dynamical correction to Dd( 31 P) À1.8 ppm, the KT2 and KT3 methods performed less accurately as compared to the B3LYP method.
3.4. The 2 J P,C NMR spin-spin coupling The 2 J P,C measured in DEP and cDEP was 5.5 Hz and 5.4 Hz, respectively. The sign of 2 J P,C was not measured; however, the negative value of 2 J P,C in phosphate was confirmed with the CCSD method. 37 The Iglo-III-calculated 2 J P,C in geometryoptimized DEP and cDEP was À7.63 Hz and À6.83 Hz, respectively. When Iglo-III was decontracted, the 2 J P,C changed only by 0.02 Hz, which illustrated the optimal contraction of the Kutzelnigg's basis (Table 4). When the FC contribution to 2 J P,C was calculated with the cc-pV5Z basis, the 2 J P,C in DEP and cDEP was À7.01 Hz and À6.31 Hz, respectively. The 2 J P,C calculated completely with the cc-pV5Z basis changed only by 0.01 Hz, which demonstrated the negligible basis set effect on the DSO, PSO and SD terms. The FC contribution calculated in DEP and cDEP with cc-pV5Z was larger by ca. 0.3 Hz than the FC with pcJ-4 (Tables S4 and S6, ESI †). The MD-averaging of 2 J P,C was calculated with Iglo-III except for the FC term, which was calculated with the cc-pV5Z basis.
The 2 J P,C in DEP and cDEP increased due to MD-averaging towards the experiment by 1.06 Hz and 0.42 Hz, respectively. The dynamic effect on 2 J P,C can be foreseen from the dependence of 2 J P,C on the T2 and T3 torsion angles (Fig. S14, ESI †) by considering the MD-calculated geometries (Fig. 2). The 2 J P,C in the snapshots fluctuated notably; however, convergence of the MD-averaged 2 J P,C with the S M smaller than 0.1 Hz was achieved within 350 snapshots ( Fig. S9 and S10). The S M deviation  represented ca. 2% of 2 J P,C . The S M for s( 31 P) represented less than 1% of the s( 31 P) magnitude (Table S8, ESI †). The numerical accuracy of the MD-averaged 31 P NMR shift was therefore comparable or even better than the accuracy of MD-averaged 2 J P,C . The effect of the explicit water solvent on 2 J P,C was calculated with 6-31+G(d) although the basis cannot guarantee the absolute accuracy of 2 J P,C ( Table 4). The 2 J P,C in geometry optimized DEP and cDEP increased due to the explicit hydration towards the experiment by 0.30 Hz and 0.58 Hz, respectively. The MD-averaged 2 J P,C in DEP and cDEP increased upon explicit hydration by 0.13 Hz and 0.34 Hz, respectively. Interestingly, a larger effect on 2 J P,C in cDEP was calculated throughout the MD simulation, which implies that the two phosphates were hydrated differently (Fig. S11, ESI †). Accessibility of the water solvent to the phosphate was assumed from the overall number of H-bonds between the two non-esterified oxygen atoms of phosphate and the hydrogen atoms of the surrounding water molecules (Fig. S15, ESI †). The average number of H-bonds for DEP and cDEP was 2.7 and 3.2, respectively. The cDEP phosphate was more exposed to solvent due to its locked ring structure while the DEP phosphate was more shielded from water due to the freely moving ethyl ends. Specific hydration of the phosphate therefore affects the magnitude of 2 J P,C . The calculated correction of 2 J P,C due to explicit hydration was 0.13 Hz and 0.34 Hz in DEP and cDEP, respectively. The MD-averaged 2 J P,C in DEP and cDEP including the hydration correction was À5.82 Hz and À5.55 Hz, respectively. The 2 J P,C increased due to geometry reoptimization with the atomic bases used in NMR calculations (Fig. 3, Tables S3 and S5, ESI †). For example, the 2 J P,C in DEP and cDEP increased due to cc-PVTZ geometry optimization by 0.45 Hz and 0.29 Hz, respectively. However, the effect on the MD-averaged 2 J P,C in cDEP was both increasing and decreasing (Fig. S12, ESI †). Consequently, the MD-averaged 2 J P,C decreased due to the cc-PVTZ geometry reoptimization only negligibly, by 0.04 Hz.
The MD-averaged 2 J P,C calculated in DEP and cDEP with the BP86 method was À5.73 Hz and À5.48 Hz, respectively.
The static value of 2 J P,C in DEP and cDEP was À6.60 Hz and À5.87 Hz, respectively. The increase of 2 J P,C due to MD-averaging with the BP86 and B3LYP methods was very similar, which allowed a reliable estimate of the dynamic effect on 2 J P,C . The average dynamical correction to 2 J P,C in DEP and cDEP was 0.97 Hz and 0.41 Hz, respectively. The corrections due to molecular dynamics and explicit hydration allowed reliable estimate of 2 J P,C with the BPW91, M06-2X and PBE0 methods based on the static NMR calculations (Table S12, ESI †). The 2 J P,C values calculated with the DFT methods ranged from À5.82 Hz to À5.42 Hz for DEP and from À5.89 Hz to À5.14 Hz for cDEP. The DFT methods gave a 2 J P,C differing from the experiment by less than 0.4 Hz. The performances of the DFT methods can be compared assuming the accuracy and the mutual difference of 2 J P,C in DEP and cDEP. Importantly, the 2 J P,C calculated in DEP was smaller than the 2 J P,C in cDEP phosphate except for 2 J P,C with the M06-2X and PBE0 methods. The B3LYP and BP86 methods provided similarly accurate 2 J P,C values in DEP and cDEP in agreement with the measured trend.
The relatively large difference between the static values of 2 J P,C in DEP and cDEP phosphates was apparently due to their notably different geometries of energy minima. The two 2 J P,C values got close in agreement with the experiment only upon dynamical averaging (Fig. 5). Hence, the very same magnitudes of 2 J P,C , which were measured in DEP and cDEP, were due to the particular dynamics of the two phosphates. The dynamic effect on 2 J P,C in DEP was twice as large as the effect on 2 J P,C in cDEP. The effect of MD-averaging on 2 J P,C can be clearly foreseen from the calculated dependence of 2 J P,C on the T2 and T3 torsions inclusive of the MD-snapshot geometries (Fig. S14, ESI †). For DEP, the 2 J P,C for the energy minimum was notably smaller than the magnitudes of 2 J P,C in the MD snapshots. For cDEP, the MD snapshots covered the areas where the 2 J P,C values were similar to the 2 J P,C values for the energy minimum. Reliable and accurate structural interpretation of the measured 2 J P,C thus inevitably must include the dynamic component. Interpretation of the measured 31 P NMR parameters employing empirical Fig. 5 The cumulative average of 2 J P,C NMR spin-spin coupling calculated with the B3LYP method in MD snapshots of (a) DEP (black square) and (b) cDEP (black circle). The 2 J P,C calculated for the global energy minimum is indicated with the red-dashed line. The measured 2 J P,C is indicated with the red line. rules and neglecting the effect of molecular dynamics might be misleading if not contradictory. The measured relative Dd( 31 P) NMR shift apparently indicated the structural non-equivalency of the two phosphates that is seemingly in contradiction with similar magnitudes of the 2 J P,C couplings although the two NMR parameters provided in fact the same structural information.

Conclusions
The benchmark for the 31 P NMR shift and the 2 J P,C NMR spin-spin coupling in phosphate was obtained by means of theoretical MD/DFT calculations considering the effects of the DFT method, atomic basis, hydration and molecular dynamics in a one-to-one relationship with the NMR measurements in two distinctively different phosphates as regards their optimal geometry and dynamic behaviour. The benchmark for the NMR parameters was aimed at the dynamically flexible DEP phosphate and the dynamically confined cDEP phosphate with notably different optimal geometries.
The calculations of the 31 P NMR shift demonstrated that the relative referencing of the 31 P NMR shift in chemically equivalent phosphates is superior to the NMR reference schemes employing H 3 PO 4 and/or a secondary PH 3 reference. The relative MD-averaged Dd( 31 P) NMR shift calculated with the B3LYP method, Iglo-III basis and PCM hydration differed from the experiment by 0.16 ppm. The calculation of 2 J P,C required a better atomic basis and inclusion of the effect due to explicit hydration that depended on the particular accessibility of the water solvent to DEP and cDEP phosphates. The MD-averaged 2 J P,C calculated in DEP and cDEP with the B3LYP method, and Iglo-III(DSO, PSO, SD)/FC(cc-pV5Z) bases inclusive of the hydration effect differed from the experiment by 0.3 Hz and 0.2 Hz, respectively.
The NMR calculations demonstrated that reliable and accurate structural interpretation of the 31 P NMR parameters in phosphate inevitably includes both the structural and dynamic components.

Conflicts of interest
There are no conflicts of interest to declare.