Interaction between Cellobiose Dehydrogenase and Lytic Polysaccharide Monooxygenase.

Lytic polysaccharide monooxygenases (LPMOs) are ubiquitous oxidoreductases, facilitating the degradation of polymeric carbohydrates in biomass. Cellobiose dehydrogenase (CDH) is a biologically relevant electron donor in this process, with the electrons resulting from cellobiose oxidation being shuttled from the CDH dehydrogenase domain to its cytochrome domain and then to the LPMO catalytic site. In this work, we investigate the interaction of four Neurospora crassa LPMOs and five CDH cytochrome domains from different species using computational methods. We used HADDOCK to perform protein-protein docking experiments on all 20 combinations and subsequently to select four complexes for extensive molecular dynamics simulations. The potential of mean force is computed for a rotation of the cytochrome domain relative to LPMO. We find that the LPMO loops are largely responsible for the preferred orientations of the cytochrome domains. This leads us to postulate a hybrid version of NcLPMO9F, with exchanged loops and predicted altered cytochrome binding preferences for this variant. Our work provides insight into the possible mechanisms of electron transfer between the two protein systems, in agreement with and complementary to previously published experimental data.

Biomass is an abundantly accessible raw material and therefore triggers continuous efforts to promote renewable energies such as biofuel. A major portion of the available biomass is composed of lignocellulose. Because of the insoluble crystalline nature of the cellulose, the glycosidic linkages are protected from hydrolytic enzymes such as cellobiohydrolases. In addition, the combination of pretreatment and enzymatic processing remains a major hurdle in the production of fermentable sugars.[1] However, the discovery that copper-dependent lytic polysaccharide monooxygenases (LPMO, CAZy auxiliary activity families AA9–AA11 and AA13–AA15) can potentiate the activity of hydrolytic enzymes by increasing the substrate accessibility opens new possibilities.[2−10] The substrate accessibility of hydrolytic enzymes is being improved by the oxidative cleavage of the polysaccharide backbone performed by LPMO.[5] The details of the underlying mechanism are still under investigation, especially the nature of the co-substrate that has been reported to be either molecular oxygen or hydrogen peroxide.[3,5,11−14] Nonetheless, studies have agreed that an electron donor is needed to reduce the active-site copper in LPMO prior to the activation of the oxygen species. Reduction of the LPMO Cu(II) center has been reported to occur over several possible routes, which include small molecular reductants of fungal and plant origin as well as photosynthetic pigments.[5,15−17] Another electron transfer system includes extracellular proteins that reduce the LPMO Cu(II) center either directly [e.g., flavocytochrome cellobiose dehydrogenase (CDH; CAZy family AA3.1)] or indirectly through redox mediators that can be recycled by glucose-methanol-choline (GMC) oxidoreductase family enzymes (CAZy subfamily AA3.2).[5,11,12,17] In the case of CDH, the electron transfer chain is starting at the dehydrogenase domain (DH) containing flavin adenine dinucleotide (FAD) where two electrons are gained by the oxidation of cellobiose into cellobionolactone.[18] In a subsequent interdomain electron transfer, one electron is being passed from the FAD to heme b of the cytochrome domain (CYT).[11,19−22] In the open state, the CYT is connected to the DH domain via a flexible linker.[11] This allows the CYT to interact largely independently with the LPMO and reduce the Cu(II).[11,17] It has been proposed that the CYT interaction site is located in the vicinity of a conserved PGP motif on the LPMO surface.[23] In such a model, the electron would follow a path from CYT heme b through the LPMO over several amino acids toward the Cu(II) center, allowing the simultaneous binding of the substrate and CYT to LPMO. Conversely, in a direct electron transfer model, the electron transfer occurs directly between CYT heme b and the LPMO Cu(II) and would not allow concurrent interaction of LPMO with its polysaccharide substrate and CYT. The latter model is supported by experimental data reporting an overlap in the interaction surface patches for the substrate and CYT on LPMO.[24] However, the amino acids involved in the interaction between LPMO and CYT were limited to a narrow surface patch around the LPMO active center that is pointing to a transient and promiscuous interaction.[24] A direct interaction between CYT heme b and LPMO copper center was also recently confirmed by small angle X-ray scattering data.[25] Here we use computer models to investigate the most likely interaction of the CYT and LPMO. We have performed protein–protein docking experiments for all combinations of a set of five related CYTs and four LPMO catalytic domains. The optimal relative orientation was further studied using molecular dynamics (MD) simulations and free energy calculations. On the basis of the results, we suggest that the interaction between the CYT and LPMO is guided by the loops that surround the LPMO active site.

Methods

Protein Structures

For the cytochrome domains, the crystal structures of Neurospora crassa CDHIIA [NcCDHIIA; Protein Data Bank (PDB) entry 4qi7], Crassicarpon hotsonii (syn. Myriococcum thermophilum) CYT (ChCYT; PDB entry 4qi3), and Phanerochaete chrysosporium CYT (PcCYT; PDB entry 1d7b) were used. In the case of NcCDHIIA, only the CYT residues (T2–C211) and heme b of the model were retained (NcCYTIIA). The amino acid sequences of N. crassa CDHIIB CYT (NcCYTIIB; UniProt entry Q7S0Y1) and Crassicarpon thermophilum (syn. Corynascus thermophilus) CDHIIB CYT (CtCYT; UniProt entry E7D6B9) were modeled onto the crystal structure of ChCYT (PDB entry 4qi6) using SWISS-MODEL.[26−28] For NcCYTIIB, residues Q1–T212 were modeled with a QMEAN4 z-score of −4.01, and for CtCYT, we used residues Q1–T203, leading to a QMEAN4 z-score of −2.94. Subsequently, the heme b cofactor was fitted into the models using the PyMOL Molecular Graphics System (version 1.7.0.0, Schrödinger, LLC). After an initial steepest descent minimization with a convergence criterion of 0.1 kJ mol–1 using the GROMOS software package for molecular simulation[29] together with the 54A7 force field,[30] MD simulations were performed as described below for 30 and 100 ns of the resulting NcCYTIIB and CtCYT models, respectively. The coordinates were clustered according to their root-mean-square deviation (RMSD) using 0.12 and 0.15 nm cutoffs to create 10 clusters for NcCYTIIB and CtCYT, respectively. The resulting clusters were used as input coordinates for the molecular docking experiments described below.

For LPMO, the crystal structures of N. crassa LPMO9C, -9D, -9F, and -9M were used: NcLPMO9C (PDB entry 4d7u), NcLPMO9D (PDB entry 4eir), NcLPMO9F (PDB entry 4qi8), and NcLPMO9M (PDB entry 4eis). To create the structure of a suggested NcLPMO9Floop variant, the three-dimensional coordinates of NcLPMO9F residues G66–P68 were replaced by those of NcLPMO9C residues F61–N77. As a further refinement, the NcLPMO9Floop structure was relaxed by steepest descent energy minimization with a convergence criterion of 0.1 kJ mol–1.

Molecular Docking

Protonation states of each molecule at pH 6.0 were calculated with PROPKA3,[31−33] and molecular docking was performed with HADDOCK 2.2.[34,35] HADDOCK distinguishes between ab initio docking and restraint docking, in which ambiguous interaction restraints can be used to guide the docking process. Furthermore, HADDOCK uses multiple stages of docking, in which an initial rigid body docking is followed by a semiflexible refinement docking stage. To limit the interaction surface of the CYT to the cofactor facing site during ab initio docking, only selected residues were made available for random ambiguous interaction restraints (Table S1) while all accessible residues of LPMO were considered for random ambiguous interaction restraints. To force contact between the two molecules, center of mass restraints were switched on, whereas the interaction between heme b and the LPMO Cu(II) center had to be satisfied during restraint docking experiments. The numbers of structures for rigid body dockings were set to 104 and 103 solutions for ab initio and restraint docking, respectively. Subsequently, 200 solutions were considered for semiflexible docking during both ab initio and restraint docking. The terms of the HADDOCK scoring function were weighted as described in Table S2. Restraint docking was repeated with three different seeds, and docking solutions were clustered according to their RMSD using a 0.15 nm cutoff. The final docking solutions were further analyzed with the GROMOS++ software package.[36]

Molecular Dynamics Simulation Settings

The optimized relative orientation of CYT and LPMO was further studied by computing the potential of mean force (PMF) of rotation around a defined axis.[37] All MD simulations were performed with the GROMOS11 biomolecular simulation package.[29] The starting coordinates were taken from the best ranked restraint docking solution for complexes with NcCYTIIA, ChCYT, NcLPMO9C, NcLPMO9F, and NcLPMO9Floop. The PyMOL Molecular Graphics System (version 1.7.0.0, Schrödinger, LLC) was used to create sets of 36 input coordinates for each pair by rotating the individual molecules in 10° steps relative to each other around the axis given by the LPMO center of mass and its Cu(II) atom. The resulting structures were parametrized using the GROMOS++ software package[36] with the GROMOS 54A7 force field.[30] Heme b was parametrized in its reduced form using the parameters described by Zou et al.[38] Partial atomic charges around the copper ion of LPMO were obtained from density functional theory (DFT) calculations on the Cu(II) center and its coordination ligands (NcLPMO9F, H1, H72, and Y157). The geometry of the LMPO fragments was derived from the NcLPMO9F crystal structure (PDB entry 4qi8); in the DFT calculation, the Cu(II) coordination was kept fixed and only atoms of ligands were optimized. The calculations were performed by means of the Gaussian package[39] using the PBE functional[40] and the def2-TZVP atomic basis set.[41] The partial atomic charges were obtained from natural population analysis[42] and are listed in Figure S1.

Subsequently, the parametrized complexes were subjected to two in vacuo steepest descent energy minimizations with a convergence criterion of 0.1 kJ mol–1. Only after the first energy minimization was the SHAKE algorithm used to constrain the bond lengths, allowing unfavorable interactions between side chains to be released during the initial minimization. The Cu(II) atom was stabilized in the LPMO active site using distance restraints with a force constant of 104 kJ mol–1 nm–2 (see Table S3 for atoms and distances), and the iron center of CYT heme b was ligated to its coordinating residues (i.e., HIS176 Nε and MET74 S) via a covalent bond to avoid any dissociation. The rotation of the individual complex units relative to each other was restrained using a dihedral angle restraint with a force constant of 0.05 kJ mol–1 deg–2, roughly allowing a 10° deviation from the target value in either direction. This restraint was defined by four points (i–j–k–l) for every complex as shown in Table S4. The latter restraints were maintained during all of the following steps. Once the proteins were solvated in a periodic rectangular box with simple point charge (SPC) water[43] maintaining a minimum solute–wall distance of 0.8 nm, the systems were energy minimized again to remove unfavorable solute–solvent contacts. Subsequently, sodium and chloride atoms were added to ensure a neutral system charge (Table S6). Thereafter, the systems were equilibrated at 50 K with initial random velocities generated from a Maxwell–Boltzmann distribution. During five discrete steps of 20 ps, the temperature of the systems was increased to 300 K. Finally, roto-translational constraints on the CYTs were initialized in two additional equilibration steps to hamper the overall rotation of the whole complex. During the production simulations, a weak coupling scheme with coupling times τT = 0.1 ps and τP = 0.5 ps and an isothermal compressibility of 4.575 × 10–4 kJ–1 mol nm3 were used to maintain a constant temperature of 300 K and a constant pressure of 1 atm.[44] The reaction field method[45] was used to treat nonbonded interactions with a cutoff radius of 1.4 nm and an ε of 61 to represent the medium outside the cutoff sphere. On the basis of the histograms shown in Figures S6–S11, additional simulations were added to ensure sufficient overlap for the computation of the PMF plot. At least 40 rotational states were simulated for at least 6 ns. Up to 2 ns was discarded for equilibration to ensure converged histograms. In addition, plain MD simulation of homology models and ab initio docking solutions were performed for 100 ns using the same settings as described above, however without the dihedral angle restraint between the two proteins. Furthermore, the force constant for the Cu(II)-coordinating distance restraints was reduced to 400 kJ mol–1 nm–2.

Analysis

The PyMOL Molecular Graphics System (version 1.7.0.0, Schrödinger, LLC) was used to visualize and analyze the docking solutions and protein complexes resulting from the MD simulations. Geometric properties between the CYT and LPMO were calculated using the GROMOS++ software package.[36] Thereafter, the resulting dihedral angle measurements were taken as input for the weighted histogram analysis method (WHAM) script implemented in the GROMACS molecular simulation suite[46] to compute the potentials of mean force (PMFs) in a periodic manner with a bin size of 1°. Error estimates were calculated using a bootstrapping procedure with 100 repeats. A Savitzky–Golay filter was applied on the PMF data with a filter window length of 69 and a polynomial of the first order.[47,48] The following python packages were used for data analys and representation: SciPy,[48] NumPy,[49] and Matplotlib.[50]

Results

Defining the Interaction Region of the CYT on LPMO

To investigate possible interaction sites on LPMOs, guided ab initio docking of 20 LPMO–CYT complexes was performed. Because the interaction surface on CYTs is given by an area surrounding their heme b cofactors,[24] the ab initio docking was guided by defining accessible residues on CYTs accordingly (Table S1). As one can see in Figure , the centers of mass (COMs) of the 1000 best ranked individual CYT positions were scattered around the LPMOs. It is only when looking at the top fifth of the 1000 best ranked complexes that COMs separate into distinct groups (Figure and Figure S2). In the case of NcLPMO9C, -9D, and -9M, a majority of the COMs of NcCYTIIA, NcCYTIIB, ChCYT, PcCYT, and CtCYT were in the vicinity of the LPMO Cu(II) centers (group 1), whereas only a small number of CYTs were located elsewhere around the LPMOs (panels A, B, and D, respectively of Figure and Figure S2). In contrast to this clear distribution, the CYT COMs around NcLPMO9F can be divided into two groups (Figure C). The first cluster was again in the vicinity of the LPMO Cu(II) center like group 1 described above, while the cluster of CYTs was in an area of NcLPMO9F that was previously proposed as a possible binding site for the CYT.[23] The ratio between the two groups was ∼1 for all CYTs, and only a small number of CYT positions could not be assigned to either group (Figure S2). In the case of the CYT positions included in group 1, the electron tunnelling pathway as computed using the gromos++ program epath[51,52] from the heme b iron to the LPMO Cu(II) center measured 1.29 ± 0.03 nm and involved a single jump through space. Conversely, the electron tunnelling pathway for CYT positions included in group 2 was significantly longer (5.9 ± 0.5 nm). Furthermore, a plain MD simulation of a NcLPMO9F–NcCYTIIA complex from group 2 showed large fluctuations in terms of contact surface area and LPMO–heme b distance over time (Figure S3), which indicates a weaker interaction in this position than with the area around the copper center. As we will show below, simulations on group 1 conformations lead to significantly smaller fluctuations in these properties (Figure S4).

Figure 1

Distribution of CYT centers of mass (COMs) around LPMOs after ab initio docking. (A) NcLPMO9C, (B) NcLPMO9D, (C) NcLPMO9F, and (D) NcLPMO9M are depicted in a yellow cartoon representation, and their Cu(II) centers are shown as spheres. The CYT COMs for NcCYTIIA, NcCYTIIB, ChCYT, PcCYT, and CtCYT (colored red, orange, blue, gray, and green, respectively) are shown as spheres around the LPMOs. The COMs of the last 800 out of 1000 ranked HADDOCK docking solutions are shown in reduced transparency.

Refining the CYT–LPMO Complex (restraint docking)

The guided ab initio docking was in agreement with the published data,[24,25] which suggested that the interaction region on LPMO is around its catalytic site. To predict possible modes of binding between the interaction partners allowing a direct electron transfer between CYT heme b and the LPMO Cu(II) center, docking experiments of the 20 CYT–LPMO combinations were performed with a distance restraint between any atom of CYT heme b and the LPMO Cu(II) with three different seeds. Subsequently, a clustering of the docking solutions revealed major differences in the rotation of CYT relative to LPMO, which could be measured by a dihedral angle between the two proteins (Figure and Table S4). The difference in rotation between the individual complexes could reach 180° [e.g., NcCYTIIA and NcLPMO9F (Figure A)].

Figure 2

(A) Two clusters of restraint docking solutions of NcCYTIIA with NcLPMO9F. CYTs of docking poses included in cluster 1 (red) are rotated ∼160° compared to CYTs included in cluster 2 (blue) around an axis given by the LPMO COM and its copper center. (B) Representation of the atoms used for the angle (Table S4). The dihedral angle i–j–k–l (green, blue, orange, and magenta spheres, respectively) is represented by yellow lines. The h–i–j–k (red, green, blue, and orange spheres, respectively) and j–k–l–m (blue, orange, magenta, and cyan spheres, respectively) dihedral angles are shown as black lines.

Additionally, the resulting rotational angles hinted at the presence of distinct groups with occasional overlap (Figure ). Furthermore, the maxima in the distribution of the rotational angle of NcLPMO9C and -9D with any of the five CYTs are well-defined and consistently found between 0° and 180°. Inversely, the maxima of the rotational angle distributions of NcLPMO9F and NcLPMO9M complexes are located between 180° and 315° apart from those of the NcCYTIIA–NcLPMO9M and the NcCYTIIB–NcLPMO9F complexes. The presence of distribution maxima suggests the presence of multiple binding modes defined by an optimal dihedral angle between the CYT and LPMO. However, docking experiments with NcCYTIIB and CtCYT exhibited a more undefined distribution of rotational angles spreading from 0° to 360°, which might suggest a less specific interaction between these two CYTs and the four tested LPMOs.

Figure 3

Dihedral angle distribution of restraint docking experiments. The i–j–k–l dihedral angle is defined as described in Table S4 (Figure B). Data for NcLPMO9F, -9C, -9D, and -9M are colored gray, red, green, and blue, respectively. The panels separate the complexes together with NcCYTIIA, NcCYTIIB, ChCYT, CtCYT, and PcCYT from top to bottom, respectively.

A heat map of the interacting LPMO residues from all of the docking poses shows that only residues around Cu(II) significantly contribute to the interactions with the CYT (Figure A–D) and that the contribution is decreasing circularly from the Cu(II) center. A similar analysis of the CYT surface revealed that only the site of the surface involving the heme b propionate groups is involved in the docking and not the site involving the vinyl groups (Figure E–I). Heme b propionate A was observed to contact the LPMO Cu(II) in 99.7% of the docking solutions as opposed to 8.5% for heme b propionate D. This is readily explained by the higher accessibility of heme b propionate A. It is important to notice that the interaction restraint was set between any atom of CYT heme b and the LPMO Cu(II) center. A conserved tyrosine located next to the heme b propionate groups of the CYT (i.e., NcCYTIIA Y99, NcCYTIIB Y105, ChCYT Y99, PcCYT Y90, and CtCYT Y97) showed a large contribution in the overall contacts, as well.

Figure 4

Heat map of the interacting residues of all of the docking poses on the CYT and LPMO surface: (A) NcLPMO9F, (B) NcLPMO9C, (C) NcLPMO9D, (D) NcLPMO9M, (E) ChCYT, (F) PcCYT, (G) NcCYTIIA, (H) NcCYTIIB, and (I) CtCYT. Residues are colored according to the frequency of their involvement in interacting with any atom in the restraint docking experiments, from red (100%) to blue (>0%). Areas colored gray were not involved in the interaction.

The Most Energetically Favorable Position Is Dependent on LPMO

Having a rather flat geometry, the interface between CYT and LPMO allows for multiple interaction modes in terms of the rotation relative to each other. To reveal the most favorable docking pose, multiple sets of MD simulations were performed in which the individual proteins were restrained to specific dihedral angles by a harmonic potential. From these, the free energy surface along the dihedral angles (i.e., the PMF) was computed via the weighted histogram analysis method.[46,53] To verify whether the rotational angle was influenced by the LPMO or CYT variant, at least 40 restraint MD simulations with four different LPMO–CYT complexes were performed. See the Supporting Information for the resulting dihedral angle distributions (Figures S6–S11) and a movie visualizing the rotation of NcCYTIIA around NcLPMO9C. Figure A shows the free energy landscape of NcLPMO9C and NcLPMO9F in complex with either NcCYTIIA or ChCYT. Although the calculated charges of NcCYTIIA and ChCYT (i.e., −3 and −15, respectively, at pH 6 with reduced heme b) differed notably, the energy minimum seemed to be influenced predominantly by the LPMO variant rather than the CYT. The complexes with NcLPMO9C (calculated charge of 0 at pH 6 with Cu2+) had a clear energy minimum at 210°, flanked by high energy barriers of 38.7 kJ mol–1 at 40° and 33.1 kJ mol–1 at 0° for ChCYT and NcCYTIIA, respectively. In contrast, the PMFs for complexes with NcLPMO9F (calculated charge of +4 at pH 6 with Cu2+) were less pronounced and shallower, with minima at 0° and 135° (ChCYT) and 206° and 15° (NcCYTIIA). Note that while the preference is inverted for ChCYT, both CYTs can be expected to bind with angles of 0° and 160°/200°. The maximum free energies were found to be 21.3 kJ mol–1 at 225.5° and 11 kJ mol–1 at 73.5° for pairs with ChCYT and NcCYTIIA, respectively. Interestingly, the most unfavorable torsional angle for the NcLPMO9F and ChCYT complex overlaps with the minimum energy region of the other three complexes and the least favorable angle of NcLPMO9C corresponds to a (local) minimum for NcLPMO9F. A closer look at a snapshot from the minimum regions of complexes with NcLPMO9C (Figure A) reveals that NcLPMO9C residues H64–S82 [loop 3 (Figure C)] form a major component of the interaction surface with CYT. Because NcLPMO9F loop 3 is shorter [P68–Y71 (Figure D)], the contribution to the interaction surface with the CYT is less pronounced and the interaction surface is reduced by 1 nm2.

Figure 5

Potentials of mean force for various complexes. Free energy as a function of the rotation of LPMO relative to CYT. Data for NcLPMO9C, -9F, and -9Floop are colored black, red and green, respectively. The solid and dashed lines represent data for complexes with NcCYTIIA and ChCYT, respectively. Shaded areas indicate bootstrap error estimates. Panel A shows the data for four wild type complexes, and panels B and C compare the results for a suggested variant to the corresponding curves of panel A.

Figure 6

Selected snapshots of NcCYTIIA (blue) in complex with (A) NcLPMO9C (red) and (B) NcLPMO9F (yellow) at dihedral angles (ζ) of 204° and 161°, respectively. NcCYTIIA is in the same orientation in panels A and B. The contact surface areas are 5 and 4 nm2 for panels A and B, respectively. Cartoon depictions of the (C) NcLPMO9C and (D) NcLPMO9F catalytic sites. The Cu(II) atom (brown sphere) is shown together with the histidine brace and the flanking tyrosine (stick representation). The loops around the catalytic site are numbered and colored from the N- to C-terminus as follows. Loop 1: black, NcLPMO9C N27–S34 and NcLPMO9F N29–S37. Loop 2: orange, NcLPMO9C C39–V49 and NcLPMO9F F43–T53. Loop 3: red, NcLPMO9C H64–S82 and NcLPMO9F P68–Y71. Loop 4: green, NcLPMO9C F108–W120 and NcLPMO9F Q106–S121. Loop 5: blue, NcLPMO9C L154–A163 and NcLPMO9F L145–A153. Loop 6: yellow, NcLPMO9C N202–P222 and NcLPMO9F N195–P209.

Furthermore, the entire orientation of the complex units was monitored. As one can see in Table S5 (Figure S5), the distance between the CYT and LPMO was stable around 4.5 nm within all of the simulations and the highest standard deviation that could be observed was 0.3 nm for the average of all of the NcCYTIIA–NcLPMO9C simulations. To ensure a proper measurement of the dihedral angle between the two proteins (i–j–k–l), the angles through atoms i, j, and k and j, k, and l were monitored (Figure B and Table S4). The angles were comparable during all of the simulations with averages of 64° and 72° for the i–j–k and j–k–l angles, respectively. However, on the basis of the calculated standard deviations for all of the simulations of the individual complexes, one can see that the i–j–k angle was fluctuating most with a maximum of 22° for the NcCYTIIA–NcLPMO9Floop complex. The values measured for the j–k–l angle on the other hand seemed to be more stable throughout the simulations with a maximum standard deviation of 10° (i.e., NcCYTIIA–NcLPMO9C and ChCYT–NcLPMO9Floop). Because the proteins were also able to rotate around the axes given by atoms i and j and atoms j and k, two additional dihedral angles (h–i–j–k and j–k–l–m) were monitored. We could observe that the LPMO showed increased mobility during all of the simulations around the i–j axis and the h–i–j–k dihedral angle was around 99° on average with a maximum standard deviation of 24°. In contrast, a maximum standard deviation of 8° around the k–l axis reflected a restricted movement of the CYT.

Moreover, the surface area and the minimum distance between the two enzymes were monitored for three simulations around the minimum energy region per complex (Figure S4). It could be observed that complexes involving NcLPMO9C and NcLPMO9Floop showed a slightly larger surface area compared to those of complexes with NcLPMO9F. Overall, the contact surface area is smaller than for the simulation of the group 2 complex (Figure S3). Also, the interquartile range is larger for simulations of NcLPMO9C and NcLPMO9Floop with NcCYTIIA (1.4 and 1.8 nm2, respectively) than for the remaining complexes. In all of the simulations, heme b propionate A was close to the LPMO, which also reflected the fact that all of the complexes were stable (Figure S4). Interestingly, in complexes with NcCYTIIA, the minimum distance from heme b propionate D to the LPMO was comparable to that of propionate A. On the other hand, in simulations with ChCYT, the minimum distance between heme b propionate D and LPMO was generally larger than for propionate A.

Swapping a Loop around the NcLPMO9F Catalytic Site Modulates the Energy Landscape

To test how the free energy profile along the rotation depends on NcLPMO9C loop 3, it was introduced into the NcLPMO9F model (NcLPMO9Floop). Performing the MD workflow described above with NcLPMO9Floop had a notable influence on the shape of the PMF plot (Figure B,C). At 31.1 kJ mol–1, the free energy barrier for the NcLPMO9Floop–NcCYTIIA complex was considerably higher than that of the NcLPMO9F–NcCYTIIA complex. Moreover, the minimum free energy region for the NcLPMO9Floop–NcCYTIIA complex was more contracted compared to that of the complex with the nonmodified NcLPMO9F and the local free energy minima at 15° and 329° disappeared. Although the overall shape of the NcLPMO9Floop–NcCYTIIA curve resembles the free energy profile of the NcLPMO9C–NcCYTIIA complex, the energy minimum is shifted by −25°. The same observations could be made for the NcLPMO9Floop–ChCYT complex; however, the changes in the energy landscape were more striking compared to those in the complex with NcCYTIIA. The energy minimum and maximum were shifted by 180° and 130°, respectively, for the NcLPMO9Floop–ChCYT complex compared to those of the nonmodified version of the complex. The energy barrier remained at 15.7 kJ mol–1. Interestingly, the energy minimum for the NcLPMO9Floop–ChCYT complex overlapped with the energy minimum of the complex with NcCYTIIA.

Conclusions

The ab initio docking experiments performed in this study confirmed that the main interaction surface on LPMO is located around the Cu(II) center. This is in agreement with recent experiments that were performed on apo-NcLPMO9C without a Cu(II) and NcLPMO9D in complex with NcCYTIIA.[24,25] Moreover, because of the relatively flat interaction surfaces on both LPMOs and CYTs, a rotation of the two units relative to each other is possible and was observed in the restraint docking experiments. This also indicates that the CYT propionate A–LPMO copper interaction is a pivotal point in the interaction with a significant contribution to binding. An analysis of the rotation of LPMOs relative to CYTs revealed that some complexes tend to group at certain rotational angles, while very diverse docking solutions were found for NcCYTIIB and CtCYT. Note, however, that for these CYTs, we used homology models, possibly adding uncertainty to the predictions of the complexes. Unfortunately, the overlap between the individual LPMO–CYT rotational groups was minimal after the docking experiments. To investigate the preferred orientations more systematically and allow for larger conformational changes in the protein structures, we performed extensive umbrella sampling MD simulations. For six different LPMO–CYT complexes, we could confirm that the most favorable positions in terms of rotation of LPMO relative to CYT were between 135° and 225°. The rather broad energy minimum region agrees with a transient and relatively unspecific mode of interaction between CDH and LPMO. Previously, the low specificity of molecular recognition has been associated with electron transfer between, e.g., cytochrome c and its partners.[54] Moreover, we could observe that a conserved tyrosine residue (NcCYTIIA Y99 and ChCYT Y99) in the vicinity of CYT heme b can make a polar contact with the heme b propionate D group. Furthermore, we could show that NcLPMO9C loop 3 is guiding the interaction with CYT by replacing this particular loop in NcLPMO9F. The loops around the LPMO catalytic site have previously been suggested to be involved in substrate specificity[55,56] but have hitherto never been related to the interaction with CDH. NcLPMO9C and -9D oxidize their substrate at the C4 position[12,17,57] and in our calculations show well-defined preferences of the CYT–LPMO orientation. On the other hand, NcLPMO9F is a C1 oxidizer and NcLPMO9M a C1/C4 oxidizer, which is attributed to differences in the LPMO loops.[17,55,56] Here we postulate that the differences in the loop architecture also affect the orientation of the CYT–NcLPMO9F complex, which becomes less defined in our calculations when loop 3 is shorter. Overall, this work provides insight into the interaction between CYT and LPMO and emphasizes the role of the LPMO loops surrounding the LPMO catalytic site in guiding this interaction.