Asymmetric dynamic coupling promotes alternative evolutionary pathways in an enzyme dimer

The importance of dynamic factors in enzyme evolution is gaining recognition. Here we study how the evolution of a new enzymatic activity exploits conformational tinkering and demonstrate that conversion of a dimeric phosphotriesterase to an arylesterase in Pseudomonas diminuta is accompanied by structural divergence between the two subunits. Deviations in loop conformations increase with promiscuity, leading to functionally distinct states, while they decrease during specialisation for the new function. We show that opposite loop movements in the two subunits are due to a dynamic coupling with the dimer interface, the importance of which is also corroborated by the co-evolution of the loop and interface residues. These results illuminate how protein dynamics promotes conformational heterogeneity in a dimeric enzyme, leading to alternative evolutionary pathways for the emergence of a new function.

www.nature.com/scientificreports/ However, we have observed considerable structural deviations between the two enzyme subunits during the PTE → AE trajectory, the magnitude of which has increased with promiscuity. We found that structural divergence was generated by the dynamic coupling between the loops and the dimer interface. Co-evolutionary analysis corroborated the functional link between the loop and interface residues. Taken together, our results highlight the role of protein dynamics via generating asymmetry between the enzyme subunits leading to alternative routes in the emergence of a new function.

Results
Changes in disorder correlate to the new activity. Protein disorder characterizes the preference of a sequence for a well-defined tertiary structure 24 . Disorder scores can be computed from the primary sequence and provide a coarse description of compactness and flexibility along the protein chain 25 , in accordance with NMR backbone dynamics 26 . Thus, disorder scores can used be to estimate the dynamic characteristics from a sequence in the absence of structural data. We applied this approach to characterize dynamics in all evolutionary intermediates during the PTE → AE conversion, and extended the analysis to those variants whose structures have not been determined. Thus, we used the IUPred method 27 to predict the disorder scores for the residues in all evolutionary intermediates 22 (R1-R22) and compared them to those in the starting variant (R0) in a pair-wise manner (Methods). First, we computed the root-mean-square deviations of the disorder scores (RMSD ID ), which characterised the magnitude of the changes in dynamics at a given point of the PTE → AE trajectory ( Fig. 2A). We observed that such coarsegrained changes in dynamics correlated to the decrease of the original phosphotriesterase activity (Pearson's correlations coefficient r = 0.87, significance p = 4.6 × 10 -7 ). Correlation with the emerging arylesterase activity was weaker, but significant (Pearson's correlations coefficient r = 0.79, significance p = 3.5 × 10 -5 ) ( Fig. 2A). These observations indicate that the functional switch PTE → AE requires systematic changes in protein motions.
Then, we aimed to probe the contributions of the mutations to the changes in dynamics. Thus, we determined the ratio of RMSD ID values, which were computed separately for mutated and not-mutated residues (Methods). As the disordered scores were computed using long flanking windows (Methods), mutations could affect the scores of those residues which were far-lying in sequence. Our results showed that mutations had a larger impact on dynamics than residues which were not changed in evolution (Fig. 2B). In addition, we observed that in the beginning of the PTE → AE conversion mutations affected dynamics more than structure, while towards the end of the evolution this trend was reversed (Fig. 2B). With emerging arylesterase activity (R1-R6), the dynamic contributions of the mutated residues were considerably larger than during specialisation for the new function (R8-R22).
In addition, we observed that overall dynamics increased with promiscuity (R1-R6, Fig. 2C), while exhibited a sharp drop at the generalist state, where original and new activities were comparable. During optimisation of the arylesterase activity (R8-R22) protein dynamics gradually increased, most likely to facilitate structural adjustments for the new function. www.nature.com/scientificreports/ Taken together, these results suggest that the emergence of the new function requires significant alterations to protein motions, which enable structural rearrangements for optimising the new activity. This is in accord with the proposed diminishing returns model of evolution 22 . Coupling between the loops and the interface. In this section, we investigated how the changes in dynamics were synchronized in the two enzyme subunits during the evolution of the arylesterase function. Thus, we computed the covariance between the disorder scores throughout the PTE → AE conversion using all evolutionary intermediates (Methods) (Figs. 3A, 3B). We observed strong dynamic couplings in the first part of the PTE → AE conversion, when the new function emerged (Fig. 3A). The dynamics of the L4 and L5 loops in particular exhibited a strong positive covariance with increasing promiscuity (R1-R8), while those of the L4 and L7 loops were anti-correlated. This suggests that L4 could serve as a clutch between L5 and L7 loops. The dimer interface exhibited a positive covariance with L4 and L5 dynamics, and a negative covariance with the L7 loop (Fig. 3A). However, we did not observe significant covariance between the loops or with the dimer interface during specialisation for the new function (R9-R22), in accordance with the decrease in dynamical changes in this part of the trajectory (Fig. 3B).
We also analysed the covariance between conformation and dynamics in variants whose structures had been determined (R0, R1, R2, R6, R8, R18, R22, Table S3). During the emergence of the new function (R0-R8) we observed a strong, positive covariance between L5 loop structure and L7 dynamics, and a negative covariance between L5 loop structure and L4 dynamics (Fig. 3C). This suggests that changes in dynamics preceded structural rearrangements of the L5 loop. In turn, structural changes of the L5 loop were concomitant with the dynamical  4), and was computed separately for the A (light grey) and B (dark grey) subunits. (C) Dynamics changes in coordination shells. The first shell comprises the active site (light blue), the second shell contains residues within 3.5 Å from the active site (dark blue), the third shell residues are within 3.5 Å from the second shell residues (light grey), and the fourth shell residues are the 3.5 Å from the third shell residues (dark grey). Changes in disorder were computed by Eq. (2). www.nature.com/scientificreports/ changes of the L7 loop. L5 structure was also strongly coupled to the dynamics of the dimeric interface (Fig. 3C). During the optimisation of the new function the structures of both L5 and L4 loops exhibited strong covariance with the interface (Fig. 3D). We observed opposite trends: while changes in dynamics of the interface induced conformational changes in the L5 loop, they decreased structural changes of the L4 loop (Fig. 3D). The L7 loop structure did not exhibit pronounced coupling during the specialisation for the new function. We performed co-evolutionary analysis on the coupling between the loop and interface residues using the GREMLIN method 28 (Methods, Fig. 4, Table S2). We found that residues of the L7 loop and those of the dimer interface exhibit significant co-evolutionary signatures (Fig. 4). In addition, we observed a co-evolution between the L4 and L5 loops, as well as between the N-and C-terminal parts of the dimer interface. These results indicate that dynamical couplings between residues of the loops and the dimer interface, which drive the loop rearrangements, have been evolutionary conserved (Fig. 4).
Taken together, co-evolutionary results corroborated the functional importance of the dynamical coupling between the loops and the interface residues. Covariance between structure and dynamics indicates that dynamical couplings induce conformational rearrangements for the emergence of a new function.
Functional promiscuity increases conformational divergence between the two subunits. Analysis of the structure-dynamics covariance indicated significant differences between the evolution www.nature.com/scientificreports/ of the two subunits ( Figure S1). In the A subunit, the L5 loop conformation exhibited positive covariance with the dynamics of the L4 and L7 loops and the dimer interface, while negative covariance was observed in the B subunit ( Figure S1). These results suggest that dynamic variations, in particular of the dimer interface induced opposite changes in the conformation of the L5 loop in the two subunits. Thus, we have compared the structures of the two subunits in all evolutionary variants with available structures (R0, R1, R2, R6, R8, R18, R22, Table S3). Indeed, we observed that the L5 loop exhibited remarkable deviations between the A and B subunits, especially in the generalist R6 and R8 states ( Figure S2). We further analysed the deviations in L5 conformation in the two subunits and focused on the distances between the loops (Fig. 5). We observed that in subunit B, the separation of the L5 and L7 loops (measured between A/G204 Cα and G273 Cα) decreased substantially during the PTE → AE trajectory, in particular around the generalist R6 and R8 states (Fig. 5, Table S4). In parallel, the L5 loop moved away from L4 in subunit B (measured between G174 Cα and S205 Cα). In subunit A, we observed the opposite trend, the L5-L4 separation decreased, while the distance between the L5 and L7 loops increased slightly with increasing promiscuity (Fig. 5, Table S4). The variations in the L5-L7 distance rationalise why stabilisation of L7 is coupled to destabilisation of L5 23 .
Taken together, we observed an asymmetry between the dynamical couplings between the loops and the dimer interface, reflecting the opposite movements of the loops in the two subunits.  Figure S3). These structures exhibit some differences in the conformation of the L7 loop, which considerably resize the active size. Although the L5 has a helical conformation in the closed state and is a loop in the open state, its separation from the L4 and L7 loops were rather similar ( Figure S3). The position of the L5 loop in the open and closed PTE structures was in between those in the two subunits of the generalist R8 state ( Figure S3). These observations indicate that different loop conformations may not be related to the optimisation of substrate binding versus the chemical step.
To obtain further insights into the different loop architectures, we assembled PTE structures from the Protein Data Bank, which represented snapshots of different evolutionary pathways in different organisms. We compared these structures to those of the R6 (PDB: 4xag) and R8 intermediates (PDB:4xay) and analysed the representative distances (L5-L7: A/G204 Cα-G273 Cα; L5-L4: G174 Cα-S205 Cα) between the loops. From the dataset of 160  Table S2. (B) The structure of the R0 starting variant (PDB:4pcp) with the significant co-evolving residues displayed (black dashed lines). The loops are coloured as L4 lime, L5 red, and L7 blue, and the dimer interface is shown in wheat. The substrate is magenta, the metal ions are displayed by cyan spheres.  (Table S5). We observed that enzymes with loop distances similar to the G chain (B subunit) R6 and R8 intermediates, also had considerable arylesterase activity (Table S5), or hydrolysed other substrates (not shown). In those enzymes structures which resemble the loop conformations in the A subunit primarily functioned as phosphotriesterases and were not observed to cleave other substrates. These results suggest that the two subunits of the enzyme dimer have different substrate preferences and catalytic activities in the generalist states. The A subunit is likely specialized on paraoxon hydrolysis, whereas the B subunit is promiscuous, with considerable arylesterase activity. Taken together, the comparative structural analysis of the subunit conformations to different PTE enzymes indicates that the deviation between the loop architecture has a functional relevance. While the loop conformations in the A subunit are compatible with paraoxon hydrolysis, those in the B subunit have increased preference for arylester hydrolysis. These results implicate that the evolution in the two subunits follows different pathways.

Discussion
The role of conformational heterogeneity in enzymatic catalysis 29 as well as in the emergence of new protein functions is being increasingly recognised 5 . Within this framework, mutations modulate the balance between different conformations, which are compatible with alternative functions 16 . This sheds light on the role of protein dynamics in protein evolution via modulating conformational ensembles by distant amino acid replacements. Changes in low-frequency structural fluctuations may also facilitate the emergence of new functions 7,30 , but their interpretation is far from being trivial 6 . Directed evolution of PTE from Pseudomonas diminuta to arylesterase 22 www.nature.com/scientificreports/ exemplifies how protein evolution engineers protein motions via freezing out fluctuations, which are unproductive for the emerging function.
Our results highlight a novel angle of how catalytic promiscuity stems from conformational heterogeneity. We have found that structural divergence between the two subunits of the Pseudomonas diminuta phosphotriesterase accompanies the evolution of the arylesterase function. Deviation in loop conformations between the two subunits increases with promiscuity. We have shown that dynamic couplings between the loops and the dimer interface, especially during the emergence of the new function, are responsible for the structural divergence (Fig. 3). These results rationalise the role of mutations, which are located at the dimer interface distant from the catalytic centre (Table S1). The functional importance of these dynamic couplings was further corroborated by co-evolutionary data (Fig. 4).
Comparative analysis using a set of PTE enzymes with structures and kinetic data available indicated that the different loop conformations represent functionally distinct states (Table S5). While the loop distances in the A subunit are compatible with phosphotriesterase hydrolysis, in the B subunit (G chain in R6 and R8 intermediates) they are similar to structures with promiscuous activities. Thus, we may conclude that the conversion to arylesterase in the B subunit precedes that of the A subunit. The considerable PTE activity of the A subunit, even at the end of the evolution, was beneficial for the reverse pathway from arylesterase to phosphotriesterase 31 . Quantifying the difference between the activities of the individual subunits requires state-of-the art QM/MM methods 32 because of the metal ions involved in catalysis 33 and the heterogeneity of the system 34,35 , which is being carried out in our laboratory using different simulation methods 32 .
Taken together, our results indicate that the conversion of PTE to arylesterase is not synchronized in the two subunits of the enzyme dimer. This is analogous to gene duplication; one copy may undergo significant sequence changes, while the other is more restricted to the original function 36 . Thus, we propose that dimerisation can promote alternative pathways in enzyme evolution via structural divergence due to asymmetric dynamic couplings.

Methods
Disorder predictions. The preference for a well-defined versus a disordered structure was computed using the IUPred program 27 . Disorder scores correlate to NMR backbone dynamical parameters 26 , thus they can be used to estimate coarse-grained dynamics at the residue level. We used the sequences of all the intermediates of the PTE → AE conversion 22,31 as an input for the IUPred long algorithm to predict the disorder scores. We have also performed predictions using the short version of the IUPred program, the Espritz NMR 37 and Dynamine 38 methods. The results were found to be consistent with each other.

Characterisation of changes in dynamics.
The following quantities were used to assess the impact of evolution on dynamics. The magnitude of the overall changes in dynamics was characterised by the root-meansquare deviation of the ID scores from the starting variant: where RMSD ID (X) is the root-mean square deviation of the disorder scores in the X evolutionary intermediate from the starting (R0) sequence. ID X i and ID R0 i are the ID scores of the ith residue in the X and R0 variant. The RMSD is computed from the scores of all residues, n is the number of residues in the protein.
The change in dynamics in given coordination shells relative to the starting variant was characterised by the sum of the differences between the ID scores of the corresponding residues.
where ID X i and ID R0 i are the ID scores of the ith residue in the X and R0 variant. m is the number of residues in the given (j) coordination shell.
The impact of mutations on dynamics in the X variant was quantified by the ratio of the RMSD ID (X) values of the mutated (m) and not mutated (nm) residues: where RMSD ID (X) is the root-mean square deviation of the disorder scores in the X evolutionary intermediate from the starting (R0) sequence. ID X i and ID R0 i are the ID scores of residue i in the X intermediate and the R0 starting variant. Residues were grouped based on whether they were mutated or not in the given round of the PTE → AE conversion. m is the number of mutated, nm is the number of not mutated residues.
Structure analysis. Structures of the evolutionary intermediates were assembled from the Protein Data Bank (Table S3). We noticed a considerable difference between the number and occupancy of the metal ions (Table S3) which are involved in the catalytic step. Structures of the evolutionary intermediates were compared to those of the starting variant (PDB:4pcp). The PTEs with open (PDB: 3a3w) and closed (PDB: 2r1n) active sites (Fig. 1) were derived from earlier studies 12 . PTE structures with different evolutionary pathways were assembled based on sequence similarity (Table S5).
The impact of mutations on structure in the X variant was quantified by the ratio of the RMSD Cα (X) values of the mutated (m) and not mutated (nm) residues: Scientific Reports | (2020) 10:18866 | https://doi.org/10.1038/s41598-020-75772-5 www.nature.com/scientificreports/ where RMSD Cα (X) is the root-mean square deviation of the Cα positions in the X evolutionary intermediate from the starting (R0) structure. P X i,Cα and P R0 i,Cα are the positions of the Cα atom of residue i in the X intermediate and the R0 starting variant. Residues were grouped based on whether they were mutated or not in the given round of the PTE → AE conversion. m is the number of mutated, nm is the number of not mutated residues.
In the analysis of loop distances, the corresponding residue numbers in the PDB structures were identified based on sequence alignment to the wild type PTE sequence.
Covariance analysis. The covariance of changes in dynamics was defined as: where ID i and ID j are the differences in disorder scores of residues i and j from the starting R0 variant.
The covariance of structure and dynamics changes were defined as: where ID i is the difference in disorder scores of residue i from the starting R0 variant. �Cα j is the difference in Cα of residue j from the starting R0 variant. Covariance matrices were calculated using the R software package (https ://www.R-proje ct.org/).

Co-evolution analysis.
Co-evolving residue pairs were identified with GREMLIN webserver 28 using the wild type PTE 22 sequence. The GREMLIN method is based on a multiple sequence alignment (MSA) by HHblits 39 with homologue sequences > 90% identity, > 75% coverage and < 25% gaps in MSA. The PTE (A0A060GYS1) coevolutionary analysis was performed using 819 sequences. Co-evolving amino acids were defined as residue pairs above a scaled score of 1 (Table S2).