Mechanistic basis of substrate–O2 coupling within a chitin-active lytic polysaccharide monooxygenase: An integrated NMR/EPR study

Significance Lytic polysaccharide monooxygenases (LPMOs) have unique catalytic centers, at which a single copper catalyzes the oxidative cleavage of a glycosidic bond. The mechanism by which LPMOs activate molecular oxygen is key to understanding copper (bio)catalysis but remains poorly understood, largely because the insoluble and heterogeneous nature of LPMO substrates precludes the use of usual laboratory techniques. Using an integrated NMR/EPR approach, we have unraveled structural and electronic details of the interactions of an LPMO from Bacillus licheniformis and β-chitin. EPR spectroscopy on uniformly isotope 15N-labeled 63Cu(II)-LPMO provided insight into substrate-driven rearrangement of the copper coordination sphere that predisposes the enzyme for O2 activation.


Supplementary Discussion -Calculations used in EPR analysis
Supplementary Discussion -DFT calculations in presence of substrate Table S1 -Structural statistics of apo-BlLPMO10A  Table S7 -EPR properties of BlLPMO10A with and without β-chitin Table S8 -EPR properties of the remote nitrogen atoms of the histidine rings of BlLPMO10A with and without β-chitin, determined experimentally and by DFT calculations Table S9 -Löwdin spin population analysis of 4-coordinate and 5-coordinate superoxide models (adding up to 200% for a total of 2 spins).
Table S10 -Calculated superoxide binding energies using cluster models (D) and (E).

Protein production and functional characterization
Production and purification of BlLPMO10A were performed as previously described (1). Copper saturation of purified protein was achieved by incubation with a 3-fold molar surplus of Cu(II)SO4 at room temperature for 30 min followed by removal of excess copper using a PD MidiTrap G-25 desalting column (GE Healthcare; Uppsala, Sweden) equilibrated with 20 mM Tris/HCl buffer pH 8.0 (2). Analysis of oxidized chito-oligosaccharides released in reactions with ball-milled shrimp shell (Pandalus borealis) α-chitin with a particle size of ~µm (Chitinor AS; Senjahopen, Norway) or ball milled squid β-chitin with a particle size of < 0.85 mm (France chitin; Orange, France) was performed using hydrophilic interaction chromatography (HILIC) as described by Loose et al. (2). The enzymatic reaction mixtures contained 1 μM of Cu(II)-

Sample preparation for NMR
Cloning, protein production using an LPMO expression cassette, and purification of the natural abundance and isotope-labeled ( 15 N and 13 C or 15 N) BlLPMO10A from Bacillus licheniformis (UniProt entry Q62YN7, residues 32-203), as well as conditions for NMR spectroscopy measurements, have been described previously (3,4). In order to obtain the apo form of BlLPMO10A, protein samples were incubated in 25 mM sodium phosphate, pH 5.5, containing 10 mM NaCl and 8 mM EDTA at 4 °C overnight, prior to exchanging the buffer to either 25 mM sodium phosphate, pH 5.5, 10 mM NaCl, or 20 mM MES, pH 5.5, using Vivaspin 6 spin-columns (5 kDa cut-off, Sartorius). NMR samples with Cu(II) were prepared and recorded in the following way. First, the buffer of a solution of 13 C-and 15 N-labeled apo-BlLPMO10A was changed to a 25 mM sodium acetate buffer, pH 5.5, 10 mM NaCl. The protein solution was then concentrated to 0.1 mM and ~450 µL, and reference spectra (1D-proton, 15 N-HSQC and aromatic 13 C-HSQC) were recorded. Then, CuSO4 was added to concentrations of 0.15 mM, and spectra were recorded with CuSO4 present.
NMR samples with Cu(I) were prepared and recorded in the following way. First, the buffer of a solution of 13 C-and 15 N-labeled apo-BlLPMO10A was changed to 20 mM MES, pH 5.5, followed by concentrating the protein solution to 0.1 mM and ~450 µL, after which reference spectra (1D-proton and 15 N-HSQC) were recorded. Then, the sample was incubated for 48 h in BBL GasPak Jar (Becton Dickinson, NJ, USA) chamber with a BBL CO2 Generator bag (Becton Dickinson, NJ, USA) to remove oxygen gas. The chamber was transferred to a Whitley A45 Anaerobic Workstation (Don Whitley Scientific Limited, UK), where the protein sample was transferred to a 5 mm NMR tube and Cu(I) was added in the form of a pellet (appr. 1 mg) of Cu(I)Cl. Finally, the tube was sealed with a rubber septum and parafilm. NMR spectra (1Dproton, 15 N-HSQC, aliphatic 13 C-HSQC, aromatic 13 C-HSQC, 15 N-edited-NOESY-HSQC and 13 C-edited-NOESY-HSQC) were recorded immediately after.

NMR spectroscopy
NMR spectra of BlLPMO10A were recorded at 25 °C on a Bruker Avance III 600 MHz or Avance III 800 MHz spectrometer, both equipped with a 5 mm Z-gradient CP-TCI (H/C/N) cryoprobe, at the NT-NMR-Center/ Norwegian NMR Platform in Trondheim, Norway, and on a Bruker Avance III 600 MHz equipped with a 5 mm Z-gradient Prodigy TCI (H/C/N) cryoprobe, at the Department of Chemistry and Biosciences, Aalborg University, Aalborg, Denmark. NMR data were processed using Bruker TopSpin version 3.5. NMR spectral analysis was performed using CARA version 1.5.5 (5). The NMR assignment of BlLPMO10A has been published elsewhere (3). For structure determination, three-dimensional 13 C-edited and 15 N-edited NOESY-HSQC spectra, as well as 1 H-1 H NOESY spectra were recorded. NOE cross-peaks were manually identified, assigned, and integrated using the NEASY program within CARA version 1.5.5 (5) 15 N-HSQC-type spectra acquired with different relaxation delays was exponentially fitted (6,7). The rotational correlation time was estimated using the T1/T2 ratio, assuming overall isotropic tumbling of the protein. Secondary structure elements were analyzed using the web-based version of the TALOS-N software spin.niddk.nih.gov/bax/nmrserver/talosn/ (8) using selected chemical shifts (N, H N ,C',C α , C β , H α and H β ).

PRE calculation
The PRE effect was calculated using the atomic coordinates for the 20-conformer ensemble of Cu(I)-BlLPMO10A (PDB ID: 6TWE) as input. Transverse (R2,para) PRE rates are described by the Solomon-Bloembergen equation (9,10): where μ 0 is the permeability of vacuum, γ is the proton gyromagnetic ratio, g is the electron gfactor, μ B is the free electron magnetic moment, s e is the paramagnetic electron spin number (s e = 1/2 for Cu(II)), ω I 2π is the Larmor frequency of the proton, and (ω) is the model-free (11) extension generalized spectral density function, defined by the following equation: where is the distance between the proton and the paramagnetic electron, 2 is the square of the generalized order parameter (see (12) for details regarding its calculation), is the overall protein rotational correlation time, and is the total correlation time defined as (τ c −1 + τ −1 ) −1 , where is the correlation time for the internal motion.
The calculated reduction in signal intensity I can then be calculated using the following equation: where R2,dia is the transverse relaxation rate in the absence of Cu(II) and t is the length of the 1 H transverse relaxation evolution during INEPT coherence transfers in 15 N-HSQC.
For BlLPMO10A, g was set to 2.13, which is the average of the calculated giso (Table 1); R2,dia was set to 12.5 s -1 , which is the average of 1/T2 ( Figure S3); τ c was set to 10.2 ns, calculated from the average T1/T2 ( Figure S3); τ i was set to 500 ps. The PRE calculation was performed using a script available at https://github.com/gcourtade/BlLPMO10A. Part of the script uses code available from https://github.com/KULL-Centre/DEERpredict.

NMR structure calculation
NOE cross-peak intensities were converted into distance restraints using the CALIBA (13) subroutine in CYANA 3.97 (13,14). Dihedral torsion angles (φ, ψ) predicted by TALOS-N (8) were included as conformational restraints, as was one disulfide bridge (Cys45-Cys56). Based on this input, the structure was calculated using the torsion angle dynamics program CYANA 3.97 (14). The structure calculation started by generating 200 conformers with random torsion angles, and the dihedral angles in each conformer were optimized using simulated annealing in 10,000 steps, to fit the restraints. The 20 conformers with the lowest CYANA target function values were energy-minimized using YASARA (15), first in vacuo, followed by using water as the explicit solvent and calculating electrostatics by applying the particle mesh Ewald method (16). In both these steps the YASARA force field (17)

Preparation of 63 Cu(II) stock
A 30% solution of H2O2 (1 mL) was very slowly added to 300 μL of conc. H2SO4 cooled in ice.
Pure 63 Cu sheets (60 mg, purchased from Sigma Aldrich) were added to this solution, which turned progressively blue. After 36 h, the solution was decanted and the undissolved copper metal was added to a freshly prepared H2O2/H2SO4 mixture. All the copper was dissolved after a further 18 h period. The two solutions were combined and bright blue crystals appeared within 48 h, which were then filtered, washed with EtOH and dried.

CW and HYSCORE EPR experiments
The apo forms of Spectral simulations were carried out using EasySpin 5.2.6 (19) integrated into MATLAB R2017a software. Simulation parameters are given in Table 1. g3 and |A3| values were determined accurately from the absorptions at low field. It was assumed that g and A tensors were axially coincident. Accurate determination of the g1, g2, |A1| and |A2| was obtained by simultaneous fitting of both X and Q band spectra. The spectra obtained upon addition of β-chitin presented a mixture of free enzyme and enzyme bound to the substrate. Spectra were collected up to 3 days after addition of chitin, but the free:bound ratio did not improve over time. The simulations of the chitin-bound form of the enzyme were performed after subtraction of the normalized BlLPMO10A or 15 N-BlLPMO10A (40% and 30% of free enzyme present in the samples, respectively) from the corresponding spectra.
The 14 N HYSCORE spectra were recorded near parallel (3090 G) and near perpendicular (3390 G) directions employing the same sequence reported above with τ = 136 ns and τ = 200 ns at 5 K or 20 K, collecting 256 data points in both dimensions. The relaxation decay was subtracted by baseline corrections (fitting by polynomials of 3-4 degrees) in both time domains, subsequently applying apodization (Hamming window) and zero-filling to 1,024 data points in both dimensions. After 2D Fourier transformation, the spectra were simulated using EasySpin (19).
The Davies ENDOR spectra were obtained using the sequence π − T − π/2 − τ − π − τ − echo with mw pulses of length tπ/2 = 32/128 ns and tπ = 64/256 ns. During time T a radio frequency (rf) pulse of length trf = 12 µs was generated by the Bruker DICE system and amplified by a 60 dB gain ENI A-500 W amplifier.

Geometry optimization
Atomic coordinates of BaAA10 were obtained from the crystal structure (PDB:5IJU, resolution 1.7 Å) from the Carbohydrate-Active enZYmes (CAZy) database) (20). Included in the truncated models were the central Cu(II) ion and 9 residues (His28, Glu68, Gln92, Ala123, Pro124, His125, Thr127, Trp187 and Phe196 (numbering starting at the first histidine is position 28; His 28 and His 125 are analogous to His32 and His121 in BlLPMO10A). Hydrogens were added to appropriate positions and the following modifications were made to decrease the computational cost of the calculations: His28 and His125 were truncated at the carbonyl carbon, which was replaced by a methyl group, Glu68 and Gln92 were truncated with methyl substitution of the Cγ, the nitrogen of the amide bond between Ala123 and Thr122 was replaced by methyl groups, Trp187 and Phe196 were truncated with methyl substitution of the Cβ, Thr127 was truncated with methyl substitution of the Cα and the methyl group of Cβ was removed. For the 'resting state' model, three water molecules were retained from the crystallographic coordinates, including the two coordinating water molecules, and the nearby 'distal' water molecule (O498, O454 and O526 respectively). To mimic substrate binding, a second model was made, where the water molecules from the crystal structure were replaced by one single equatorially coordinating water molecule. In order to account for the structural constraints imposed by the protein, multiple atoms were kept frozen throughout the optimization; these atoms are denoted by asterisks in Fig

EPR Property Calculations
EPR property calculations were performed on the optimized geometries using the ORCA 4.1.0 program at the DFT level of theory. The cartesian reference system was oriented as such that the NH2-Cu(II)-O axis was aligned with the y-axis and the N-Cu(II)-N axis was oriented along x.
The integration grid size was kept large (AngularGrid = 7 for all atoms and IntAcc = 7 for the Cu(II) ion) to ensure that the core density was correctly described. Solvation effects were accounted for in the property calculations by implementing the conductor-like polarized continuum model (CPCM) with a dielectric constant of 80.0 and a refractive index of 1.33 (water). The hyperfine coupling calculations included the Fermi-contact, spin dipolar and spin orbit contributions.
A variety of basis set combinations and functionals were employed to ensure extracted trends were not functional specific or basis set limited ( Table S3). The fraction of Hartree-Fock exchange was also altered in the hybrid functional, following several studies showing improvements in the EPR property calculations (22). Calculations were performed using nonrelativistic and scalar relativistic (zeroth order regular approximation, ZORA) approximations.
The overall accuracy of the spin-Hamiltonian parameters determined by DFT showed to be limited (Tables S4-S6). However, the trends predicted between Models A and B showed to be consistent. All the calculations show a change in both g and Cu(II) hyperfine tensors from rhombic to axial, going from the 'resting state' to the '4-coordinate' model. The best agreement of the experimentally derived spin-Hamiltonian parameters was achieved using scheme 2 ( Table   S3). This scheme employed the hybrid functional, B3LYP, with an adjusted fraction of Hartree-Fock exchange (38%). This scheme also utilises the IGLO-III basis set which has additional flexibility in the core region, making it more suited for EPR properties calculations with respect to the Def2-TZVP basis set. coordinate Cu II model (F) was generated to represent the geometry following superoxide release.

Superoxide model geometry optimizations
A single superoxide molecule, model (G), was also produced as the remaining product following superoxide dissociation. All amino acid residues were retained in keeping with the previous models. Geometry optimizations for models (D) and (E) were performed on the triplet potential energy surface. The individual models following superoxide dissociation (B), (F) and (G) were all optimized on the doublet potential energy surface. Both oxygen atoms on the superoxide molecules were treated using Ahlrich's Def2-TZVP basis set. All geometry optimizations were performed using the same basis set and functional schemes as outlined in the Geometry optimization section.
The atomic coordinates of all models used in the calculations are provided in the Appendix.

Superoxide bond strength calculations
To evaluate the superoxide binding energy, single point calculations were performed (as In addition, differences in the Cu(II)-superoxide bonding are evident when assessing the Löwdin spin population analyses of the two superoxide models (Table S9). These values were obtained from the same single point frequency calculation described above. The spin population on the copper ion is shown to decrease by 12.5% upon changing from 5 to 4 coordinate (mimicking substrate binding). This decrease in spin population on the metal is in accord with an analogous increase of the spin population located on the two oxygen atoms. Therefore, this analysis shows that the decrease in coordination number from 5 to 4 leads to a large increase in covalency of the Cu(II)-superoxide bond. A reduction in spin population of d(z 2 ) character (ca. 1.3%) is also seen following a reduction in coordination number.

Supplementary Discussion -Calculations used in EPR analysis
The Cu d based molecular orbitals are written as: The 2 − 2orbital is the SOMO; the represents the metal d-orbital contribution to the molecular orbital, while and are the coefficient for the 2 − 2 and 2 orbitals in the ground state orbital (GS), with 2 + 2 = 1.
The g values can be expressed as: Cu represent the one-electron quasi-atomic copper spin-orbit coupling constant (usually taken as -830 cm -1 ), and the ΔE values are excitation energies of the ligand field transitions.
Similarly, the Cu hyperfine coupling can be written as: d =gegCuμeμCu is the quasi atomic parameter usually taken as 1180 MHz, the term -PdK (in blue) represents the isotropic Fermi contact ( Fermi ), indicated in blue; the Spin-Dipolar ( SD ) contribution is in green and the Spin-Orbit contribution is indicated in red ( SO ).
The different Cu hyperfine contributions, the % 2 in the SOMO and the GS 2 were determined as follows: 1) The ratio between / (and therefore the % 2 in the ground state) was estimated from the rhombicity parameter ( ) and the g values equations following Gewirth et al.:(23) Assuming 2 ≈ 2 and Δ ⟶ 2 − 2 ≈ Δ ⟶ 2 − 2 .
2) SO was obtained from the experimental ∆ and the and values obtained in 1) 3) SO was subtracted from Total to get Fermi + SD 4) The x, y and z components of Fermi + SD were averaged to obtain Fermi 5) Fermi + SO was subtracted from Total to get SD  (Figure S7-A). The optimized geometry showed minimal distortion from the crystallographic coordinates (Figure S8), showing that the model faithfully represented the Cu(II) active site.
The spin-Hamiltonian parameters obtained from the EPR of BlLPMO10A after addition of substrate suggested that the copper coordination sphere is near axial and, in accordance with previous studies, it is reasonable to think that water molecule in the nominal axial position is displaced by binding of substrate. Hence, a second model (Model B) was produced with a square planar-like symmetry about the copper center, achieved by the removal of the pseudo-axial water molecule from the optimized geometry structure of Model A (Figure S7-B and S8-B).
Additionally, a third model was constructed based on Model B, but also including a short chitin oligosaccharide (NAG2) bound to the active site ( Figure S9).    Figure S5. The numbers in brackets represent the error on the measurement estimated from the quality of simulated fits. The Euler angles define the zy'z'' rotations with respect to the g matrix.        Table S10. Calculated energies of superoxide cluster models (D) and (E) and the resulting geometries following superoxide dissociation (B), (F) and (G). Cu(II)-superoxide binding energy calculated by the difference in electronic energy.  intensity (see Figure 1E) and within a 12 Å radius from the expected Cu(II) coordination site are colored pink, whereas residues with less than 30% remaining intensity and further than 12 Å from the Cu(II) site (Ala160-Arg162) are colored blue. The side-chains of His32 and His121 are shown as sticks.       F. Neese, Sum-over-states based multireference ab initio calculation of EPR spin Hamiltonian parameters for transition metal complexes. A case study. Magn. Reson. Chem. 42, S187-S198 (2004).