Two-step mechanism of J-domain action in driving Hsp70 function

J-domain proteins (JDPs), obligatory Hsp70 cochaperones, play critical roles in protein homeostasis. They promote key allosteric transitions that stabilize Hsp70 interaction with substrate polypeptides upon hydrolysis of its bound ATP. Although a recent crystal structure revealed the physical mode of interaction between a J-domain and an Hsp70, the structural and dynamic consequences of J-domain action once bound and how Hsp70s discriminate among its multiple JDP partners remain enigmatic. We combined free energy simulations, biochemical assays and evolutionary analyses to address these issues. Our results indicate that the invariant aspartate of the J-domain perturbs a conserved intramolecular Hsp70 network of contacts that crosses domains. This perturbation leads to destabilization of the domain-domain interface—thereby promoting the allosteric transition that triggers ATP hydrolysis. While this mechanistic step is driven by conserved residues, evolutionarily variable residues are key to initial JDP/Hsp70 recognition—via electrostatic interactions between oppositely charged surfaces. We speculate that these variable residues allow an Hsp70 to discriminate amongst JDP partners, as many of them have coevolved. Together, our data points to a two-step mode of J-domain action, a recognition stage followed by a mechanistic stage.


Introduction
By transiently binding many different polypeptide substrates, Hsp70 chaperones assist in diverse cellular processes, from de novo protein folding to protein trafficking to disassembly of protein complexes.J-domain protein (JDP) cochaperones are crucial players in the cycle of interaction of Hsp70 with substrate.They transiently interact with ATP-bound Hsp70 via their defining J-domain.This interaction, in coordination with substrate binding, facilitates stimulation of Hsp70's ATPase activity, which in turn drives the large scale conformational changes that stabilize substrate interaction [1][2][3][4].
The structure of both J-domains and Hsp70s are conserved.J-domains are composed of four α-helices.Helices II and III form a finger like structure.The helix II/III connecting loop includes a conserved histidine, proline, aspartate tripeptide (HPD).These invariant residues are critical for J-domain stimulation of Hsp70's ATPase activity.Hsp70s are more complex in structure, containing two large domains-a nucleotide binding domain (NBD) and a substrate binding domain (SBD)-connected by a "linker" segment.When ATP is bound, the domains are docked such that the substrate binding site in the β subdomain of the SBD (SBDβ) is easily accessible [5,6].Upon ATP hydrolysis, the domains disengage, and the substrate is "trapped" as the α subdomain of the SBD (SBDα) closes over the SBDβ substrate binding site [7][8][9].
The heart of the allosteric control mechanism behind these Hsp70 conformational transitions is an interaction network at the interdomain interface [10][11][12].A critical segment of this interface is formed by interactions at the base of the NBD with the linker and SBDβ (called NBD/SBDβ,linker throughout) [5,6,13].The linker plays an important role, as its interaction with the NBD is key to stimulation of ATPase activity [14][15][16][17].A transient allosterically active intermediate poised for hydrolysis of ATP, often referred to as the "allosterically active state", has been observed in the presence of a peptide substrate-NBD/SBD contacts are largely absent, but the linker remains bound to the NBD [10,18].
It has been clear for years that J-domain interaction is essential for stimulation of Hsp70's ATPase activity [19].Particularly relevant to this report, early work on Escherichia coli JDP DnaJ and Hsp70 DnaK established an arginine at the Hsp70 NBD/SBDβ,linker interface, R167 of the NBD, as an important for J-domain function [14,20,21].More recently a crystal structure of the J-domain of DnaJ in complex with ATP-bound DnaK was obtained [22].This normally transient interaction was trapped by covalently linking the J-domain and Hsp70 in cis, and by crosslinking Hsp70's NBD and SBD to stabilize intramolecular interactions between domains.Consistent with the results of a wealth of mutational and biochemical analyses, as well as earlier molecular dynamics simulations [23], the J-domain is bound at the NBD/SBDβ, linker interface.Residues of helices II and III and the connecting HPD loop form polar and hydrophobic interactions with the NBD and SBDβ.
Despite the accumulated information about J-domain/Hsp70 interactions, important questions remain.Although the end result of ATP hydrolysis is destabilization of Hsp70 domaindomain interactions, the initial structural consequences of J-domain binding are not known.Furthermore, there is no understanding of what determines the specificity of J-domain/Hsp70 interactions and how Hsp70s discriminate amongst their multiple JDP partners.To address these issues, we investigated two JDP/Hsp70 pairs: DnaJ/DnaK of E. coli and the mitochondrial Hsc20/Ssq1 of Saccharomyces cerevisiae [24,25], which specializes in the biogenesis of proteins containing iron-sulfur clusters.
Results of our analyses point to two previously unappreciated facets of J-domain/Hsp70 interactions, leading us to propose a two-step model of J-domain action-recognition, followed by mechanistic action.Molecular simulations point to an active role of the invariant aspartate of the J-domain HPD in the mechanistic step-perturbing the intramolecular interaction network formed by the previously identified, key NBD arginine.This perturbation promotes allostery-related conformational changes of Hsp70.Evolutionary analyses indicate that J-domain/ Hsp70 recognition, is dependent on variable, coevolving residues, even though the general mode of interaction of J-domains with Hsp70s is the same.This variability provides a plausible explanation for how Hsp70s are able to discriminate amongst their multiple JDP partners, thus enabling evolution of complex JDP/Hsp70 interaction networks.

Structural model of the Hsc20-Ssq1 complex
To obtain a structural model of the JDP-Hsp70 complex formed by Hsc20 and Ssq1, we combined protein-protein docking with all-atom molecular dynamics (MD) simulations (see Methods and S1 Fig for details).This approach revealed a dominant bound state that accounted for 56% of the bound population (Fig 1A,S2 Fig).No other bound state accounted for more than 6% of the ensemble; these minor states rapidly interconverted and had a clear tendency to converge to the most populated bound state (S2 Fig) .The binding mode in this dominant state closely resembles the recently published X-ray structure of DnaJ's J-domain in complex with DnaK (DnaJ JD -DnaK) [22].Helices II and III are oriented very similarly in the two complexes-with a backbone root mean square deviation (RMSD) of 4.17 Å (Fig 1B).More specifically, the J-domains bind at the NBD/SBDβ,linker interface such that helix II predominantly interacts with the β-sheet region of NBD subdomain IIa and helix III with SBDβ (Fig 1).The HPD residues are positioned nearly identically in the two complexes, towards R167/R207 of NBD subdomain Ia in DnaK/Ssq1 (referred to as R NBD throughout), the arginine that previous analyses of DnaK established as important for allosteric transitions [14,20,21].The predicted model of the Hsc20-Ssq1 complex was additionally validated using biochemical experiments, which are described later in the text.

D HPD interferes in Hsp70 interdomain interactions by perturbing critical R NBD contacts
To address the question of the mechanism of J-domain action, we employed MD simulations that, unlike most other methods, allow studying of transient residue-residue interactions [26].We first analyzed Hsp70 alone, starting with the Ssq1 structural model (see Methods for details)-focusing on intramolecular interactions involving the R NBD because of its demonstrated importance in NBD/SBD communication [14,20,21] and its proximity to the HPD in the J-domain-Hsp70 complexes.Our initial analysis of MD trajectories revealed R NBD to be capable of forming two ion pairs, one with D517 of the SBDβ (referred to as D SBD throughout) and one with D429 of the interdomain linker (referred to as D Lk throughout).To analyze these two contacts of R NBD more closely, we constructed the free energy map governing their formation (Fig 2A ), by applying metadynamics, an enhanced sampling method suitable for simulating complex conformational changes [27,28].As relevant coordinates, we used the distances between the centers of mass of the guanidine group of R NBD and carboxylic groups of its interaction partners D SBD and D Lk .Using the distance threshold for an ion pair of 0.5 nm, we found the probability of the two contacts occurring at the same time to be 82% (Fig 2A , S1 Table ).In the DnaK crystal structure [5] the homologous D SBD -R NBD pair is present, but the homologous D Lk -R NBD pair is not.To address this apparent discrepancy, we performed metadynamics simulations of the DnaK crystal structure.Analysis of the free energy map revealed that DnaK R NBD interacts simultaneously with D SBD and D Lk 83% of the time.Taken together, our metadynamics simulations indicate that the R NBD double interaction can occur frequently in both Ssq1 and DnaK, despite the fact that it was not captured in the static DnaK crystal structure.
To assess the effect of J-domain binding to Hsp70 on this conserved arginine-centered interaction network, we carried out the same set of calculations for Ssq1 and DnaK in complex with Hsc20 and DnaJ JD , respectively (Fig 2A, S1 Table ).For both Hsc20-Ssq1 and DnaJ JD -DnaK, the fraction of the population of Hsp70 molecules in which R NBD contacts both D SBD and D LK was reduced compared to that of Hsp70 alone: Ssq1, from 82 to 62%; DnaK, from 83 to 51%.We hypothesized that D HPD competes with D SBD and D Lk for interactions with R NBD and therefore is responsible for the observed partial loss of its intramolecular contacts.We also used metadynamics trajectories to calculate all R NBD contacts both in the presence and in the absence of the J-domain (S3 Fig) ; beside the aspartates discussed above, only interaction with the Ala residues next to D SBD , was reduced in the presence of the J-domain.To determine the likelihood of such an exchange of interaction partners, we performed metadynamics simulations to compute the free energy profiles for transition of R NBD from intramolecular interaction with D SBD and D Lk to interaction with D HPD (Fig 2B).In these simulations, the HPD loop was kept near NBD subdomain Ia in the orientation predisposing it for the interaction with R NBD (see Fig 1A).The well-defined free energy minima in the resulting profiles indicate that R NBD shows preference (1.5 and 4 kcal/mol for Ssq1 and DnaK, respectively) for interaction with D HPD , rather than simultaneous interaction with the aspartates on the SBD and the linker, D SBD and D Lk .The observed differences in the depth of the free energy minima between Ssq1 and DnaK could be attributed to differences in the amino acid sequences or to the fact that Hsc20 was modeled as a whole instead of only its J-domain.Finally, an example of spontaneous transition of R NBD to the intermolecular contact with D HPD , captured in our unbiased equilibrium MD simulation of the Hsc20-Ssq1 complex, is shown in S1 Movie.Taken together, our data suggest that D HPD perturbs the R NBD -centered interaction network at the NBD/ SBDβ,linker interface.

J-domain binding initiates SBDβ disengagement from the NBD
We next asked how perturbation of intramolecular interactions of R NBD with D SBD and D Lk affects the NBD/SBDβ,linker interface.First, we considered the effect of substitution of R NBD , complex to illustrate the orientation of the J-domain relative to the NBD subdomains (IA, IB, IIA, IIB s) which are marked, only the J-domain of Hsc20 is shown for clarity.(top right) Close-up showing J-domain interaction with an orientation similar to that at left.(bottom left) X-ray structure of the J-domain of DnaJ (DnaJ JD ) in complex with DnaK (DnaJ JD -DnaK) (PDB ID 5NRO).(B) Structural alignment of the predicted Hsc20-Ssq1 complex and DnaJ JD -DnaK crystal structure (the backbone RMSD of helices II and III of the J-domains after aligning NBDs is 4.17 Å).The Jdomain, NBD, linker and SBD of aligned proteins are represented in cyan/magenta, grey/pink, yellow/green and light/dark orange for Hsc20-Ssq1 and DnaJ JD -DnaK, respectively.https://doi.org/10.1371/journal.pcbi.1007913.g001such that interactions with aspartates are prohibited.An R NBD ->D DnaK substitution variant has been shown to be defective in interdomain communication [14].Its intrinsic ATPase activity is somewhat higher than that of wild-type (WT) DnaK and not effectively stimulated in the simultaneous presence of its JDP partner DnaJ and protein substrate.We therefore asked whether substitution of R NBD ->D in Ssq1 had a similar effect.We found the basal ATPase activity of Ssq1 R NBD ->D to be 3-fold higher than that of WT Ssq1, and not efficiently stimulated in the presence of Hsc20 and protein substrate, reaching only 30% of the activity observed for WT Ssq1 (S4 Fig).Thus, substitution of the conserved R NBD has significant effects on interdomain communication in both DnaK and Ssq1.
Because of these effects on interdomain communication, we used MD simulations to probe how weakening of the cross-domain interactions of R NBD (i.e. with D SBD and D Lk ) affects the compactness of the NBD/SBDβ,linker interface.Umbrella sampling was used because it allows testing of conformations not accessible by conventional MD simulations.For these simulations a center of mass distance between SBDβ and its NBD interacting region (SBDβ-NBD distance) was used as a relevant coordinate.To accelerate the free energy convergence and to allow exhaustive sampling of interdomain contacts these simulations were performed using Ssq1 and DnaK with a truncated SBDα subdomains.
Free energy profiles were computed (S5 Fig) , and then converted into probability distributions (Fig 3A).These distributions show that Ssq1 and DnaK both adopt conformations with longer SBDβ-NBD distances in the presence, than in the absence of a J-domain, indicative of a transition to a less compact interdomain interface.Furthermore, the probability distributions of the R NBD ->D variants (Fig 3A ) resemble those of Hsp70s in the presence of the J-domain in that that the mean distance values are very similar for both distributions, supporting our hypothesis that disruption of the arginine centered interaction network affects compactness of the interdomain interface.
To assess what structural changes took place within the NBD/SBDβ, linker interface we used umbrella sampling trajectories to determine the number of interdomain contacts (Fig 3B).We analyzed interactions between the NBD and SBDβ and between the NBD and the linker separately, because of their opposing effects on ATPase activity (i.e.ATPase activity requires linker association and SBD dissociation).In all our analyses, for both Ssq1 and DnaK, the number of the NBD/linker contacts remains virtually the same regardless of the interdomain distance, while the number of NBD/SBDβ contacts decreases steadily with increasing distance (Fig 3B).Taken together, these data indicate that in the presence of either the Jdomain or the R NBD ->D substitution (i.e., when the interdomain distance increases), both (a) Free energy landscape governing the formation of Ssq1/DnaK R NBD interdomain contacts with D Lk and D SBD (top) Schematic depicting of the three Hsp70 residues of the R NBD interaction network; R NBD of subdomain Ia (blue dot) can interact with D LK of the linker (yellow dot) and D SBD of the SBDβ (brown dot).The four possible interaction states of R NBD are illustrated at right with the color of the rectangles serving as a key for the free energy landscapes at bottom.NBD subdomains Ia (light grey) and IIa (dark grey), SBD subdomains α and β (brown), linker (yellow).R NBD = R207/R167 in Ssq1/DnaK; D Lk = D429/D393 in Ssq1/DnaK; D SBD = D517/D481 in Ssq1/DnaK.(bottom) Free energy landscapes governing the formation of these interaction states for Ssq1 and DnaK alone (left panels) and for Hsc20-Ssq1 and DnaJ JD -DnaK complexes (right panels) with distances between RNBD-D Lk and between R NBD -D SBD used as relevant coordinates.Integrated populations of these interaction states are shown by pie charts (see S1 Table for data  DnaK and Ssq1 adopt conformations with markedly reduced number of NBD/SBDβ contacts.In particular, upon J-domain binding, the populations of DnaK and Ssq1 that retain at least 70% of interdomain contacts are reduced from 88% to 50% and from 91% to 45%, respectively.A similar reduction is observed for the R NBD ->D substitution in the absence of J-domain, indicating that perturbation of the arginine-centered interaction network promotes partial disengagement of SBDβ.These results also indicate that J-domain binding perturbs the R NBD intramolecular interactions, which in turn leads to partial disengagement of SBDβ, while the interdomain linker remains tightly associated with the NBD.Such an arrangement favors ATP hydrolysis, which triggers further conformational changes [10,11,14,15].

Conservation of residues involved in the mechanistic step of J-domain action
The R NBD -centered interaction network is the same in DnaK and Ssq1, in that R NBD interacts with homologous aspartates in the linker and SBDβ.To investigate the prevalence of these residues, we analyzed Hsp70 sequences from a large number of bacterial and eukaryotic species.R NBD and its interaction partner D Lk is invariant in our data set (Table 1).The conservation of its SBDβ interaction partner D SBD is somewhat more complex.Aspartate is only present at this position in DnaKs of bacteria closely related to E. coli and in mitochondrial Hsp70s of fungi; asparagine is present at this position in most Hsp70s (Table 1, S6 Fig) .That asparagine can also form a hydrogen bond with arginine [29], thereby participating in the interaction network, supports the idea that the arginine-centered network is conserved across Hsp70s.Thus, it is likely the mechanism of action of J-domains, that is perturbation of this network, is conserved as well.

Importance of electrostatic interactions at the Hsc20/Ssq1 binding interface
Having found similarities of the R NBD network and the overall mode of interaction of DnaJ JD -DnaK and Hsc20-Ssq1, we decided to look more closely at the J-domain/Hsp70 binding interfaces.DnaJ JD -DnaK has been studied extensively [22]; both polar and hydrophobic contacts were identified based on the DnaJ JD -DnaK crystal structure and the functional importance of key residues was verified experimentally.However, no information is available for Hsc20-Ssq1.We therefore expanded our computational analysis of the Hsc20-Ssq1 complex.We used MD simulation trajectories to search for all specific (polar and hydrophobic) residue-residue contacts, and computed the relative energetic contributions of individual residues to complex stability (S7 Fig, S8 Fig) .The binding interface is formed by 17 residues of the J-domain contacting 18 residues of Ssq1 at a cutoff distance of 0.5 nm.The interface, which occupies 6.7 ± 0.2 nm 2 , is mostly hydrophilic, with scattered hydrophobic interactions contributing 12.5% of the surface (0.8 ± 0.1 nm 2 ).The predominant contributors are eight charged residues (Fig 4A).Three positively charged residues of Hsc20 helix II (R37 JD , K38 JD , R41 JD ) that form a  network of electrostatic interactions with four negatively charged residues of the Ssq1 NBD (D246, E248, D249, E253) and K70 JD of helix III that promotes complex formation, though its binding partner(s) was not identified in our analyses.We carried out two types of biochemical experiments to verify the functional importance of key residues that based on our computational analysis form the Hsc20/Ssq1 binding interface.First, we asked if alanine substitution of these Ssq1 residues disturbs the Hsc20/Ssq1 interaction.To enable detection of the inherently transient J-domain/Hsp70 interaction, we used Hsc20 having a site-specific p-benzoyl-l-phenylalanine (Bpa) crosslinker [30] replacing two residues, one upstream (Q46) and one downstream (M51) of the HPD (Fig 4B, S9 Fig) .The efficiency of Hsc20 crosslinking to the Ssq1 variants was strongly reduced compared to that of WT Ssq1.As a control for crosslinking specificity we replaced a residue outside the binding interface with Bpa (S34Bpa); no crosslink band(s) was detected (S10 Fig) .Second, we measured Hsc20's ability to stimulate the ATPase activity of Ssq1.We tested four Hsc20 and four Ssq1 variants, each having a different single alanine substitution.All variants had reduced activity (Fig 4C).We conclude that the key interfacial residues of Hsc20 and Ssq1 identified by the MD simulations are indeed functionally important for the Hsc20/Ssq1 interaction.

Hsc20 binding to Ssq1 is driven by long-range electrostatic interactions
The functional importance of charged residues in Hsc20/Ssq1 interaction led us to consider whether the initial encounter of the two proteins is electrostatically driven-an idea furthered by the presence of a patch of uniformly positive electrostatic potential around helix II of the Jdomain and a complementary patch of negative potential around subdomain IIa of the Ssq1 NBD, the site of J-domain interaction (S11 Fig) .MD simulations in which Hsc20 and Ssq1 were initially separated with their charged patches facing each other, revealed that Hsc20 spontaneously interacted with Ssq1 within tens of ns and remained stably bound for the rest of the stimulation (S2 Movie).Notably, the spontaneously formed complexes are very similar to the dominant bound mode identified in our systematic search (S11 Fig).
To examine the driving forces behind spontaneous Hsc20/Ssq1 interaction, we computed the underlying free energy profile using umbrella sampling.The resulting profile (S11 Fig, shows a pronounced bound-state energy minimum (at distances < 2 nm), indicative of sitespecific binding.The non-zero slope of the WT free energy profile up to ~3.5 nm demonstrates that Hsc20 is attracted to Ssq1 by long-range electrostatic driving forces.Consistent with the above analysis of the binding interface and biochemical data (Fig 4), the free energy minimum, and hence binding affinity, almost completely disappears in the case of the Hsc20 variant having helix II alanine substitutions R37A,R41A (S11 Fig).

Residues forming the JDP/Hsp70 interface are variable and coevolving
The two divergent J-domain/Hsp70 partnerships analyzed here show remarkable structural similarity in their overall mode of interaction-J-domain helix II with the β-sheet of NBD subdomain IIa and J-domain helix III with SBDβ (Fig 1).However, when we aligned amino acid sequences of these regions it became apparent that the residues forming the two J-domain/ Hsp70 binding interfaces are not evolutionary conserved (Fig 5A).This is particularly evident for the J-domains-only 12 of the 17 positions forming the Hsc20/Ssq1 interface are shared with the DnaJ JD /DnaK interface, and of these 12, only 3 are occupied by the same amino acids in the two interfaces.A broader comparison of Hsc20, Ssq1, DnaJ, and DnaK with their orthologs from a wide range of species belonging to fungi, animals and plants (S12 Fig, S13 Fig) revealed that the amino acids occupying the J-domain/Hsp70 interfacial positions are indeed variable.The J-domain side of these interfaces is particularly so.Amino acid sequence variability is evident for J-domains of Hsc20 and DnaJ among relatively closely related species, Ascomycota and Bacteria, respectively, suggesting that even closely related JDP/Hsp70 pairs engage different sets of residues to form the binding interface (Fig 5B).One trend is evident however.Despite the overall sequence variability, the J-domain side of the interface is enriched in positively charged amino acids (Fig 5C )-consistent with our results that complementary electrostatic interactions are crucial for stability of the Hsc20-Ssq1 complex.The Hsp70 side of   the binding interface is less variable than the J-domain side; 8 out of the 16 positions shared between the Ssq1/DnaK interfaces are occupied by identical amino acids.
The data discussed in the previous paragraph presents an apparent conundrum-sequence variability of J-domain/Hsp70 interfacial residues, yet maintenance of the functional interaction of J-domains with Hsp70 partners.Coevolution of residues across the binding interface, such that amino acid substitutions on one side are compensated by complementing substitutions on the other [31], could be an explanation.
Therefore, we tested whether J-domain/Hsp70 interfacial residues are coevolving.We analyzed the Ascomycota and bacterial data sets discussed above (Fig 5B ), applying the probabilistic Coev model [32] of sequence evolution that allowed us to discriminate between coevolving and independently evolving positions, while incorporating the underlying phylogenetic relationships among analyzed sequences (S14 Fig, S15 Fig) .We found that some of the variable positions are indeed coevolving-16 for the Hsc20/Hsp70 interface; 20 for DnaJ/DnaK (Fig 5B).As a control, we tested the same datasets using a model of 'independent' sequence evolution.The Coev model fits our data significantly better than the independent model (S2

Discussion
Based on the results presented we propose a mechanism of J-domain function that reconciles two apparently contradictory features: an invariant dependence on the HPD signature motif, yet an ability to recognize a specific Hsp70 partner.Our data supports a model in which the D HPD of the J-domain plays an active role in perturbing the intramolecular interactions of the Hsp70 R NBD with conserved residues of the linker and SBDβ, destabilizing NBD/SBD interactions and thus promoting allosteric transition leading to ATP hydrolysis.On the other hand, evolutionarily variable J-domain residues are responsible for the specificity of exclusive JDP/ Hsp70 interactions, as well as fine-tuning prioritization of JDP/Hsp70 interactions in the cases where an Hsp70 has multiple JDP partners, each of them guiding Hsp70 to different cellular functions.

Mechanistic step
Based on detection of an energetically favorable interaction between invariant D HPD and invariant R NBD , we propose that the DHPD-R NBD interaction displaces those of R NBD with residues of the interdomain linker (D Lk ) and SBDβ subdomain (D SBD ) of Hsp70.This displacement promotes destabilization of the interdomain interface, such that the SBDβ partially disengages from the NBD.These results are consistent with experimental data showing less protection against hydrogen exchange of this region when residues in the NBD/SBDβ interface are mutated [33] It is important to note that the interdomain linker remains stably bound to the NBD upon J-domain induced disengagement of the SBDβ.This structural arrangement has been termed the "allosterically active state" [10] because linker engagement with, but SBDβ disengagement from, the NBD is necessary for ATPase activation that triggers further Hsp70 conformational changes required for stabilization of substrate binding [11, 14-17, 34, 35].The data presented here also allows a more unified picture of how Hsp70 reaches the allosterically active state to effectively trap substrate.It was previously shown that the allosterically active state is populated in the presence of ATP and substrate [10,11,36].Substrate binding induces an allosteric signal, which is propagated through a cascade of interactions from the SBDβ substrate binding site to the NBD/SBDβ,linker interface.Our data indicates that the allosterically active state it is also populated in the presence of ATP and J-domain.In this case, the signal is propagated by perturbation of the R NBD centered interaction network.These perturbations could also modify the interaction mode of ATP bound to Hsp70, thus affecting its rate of hydrolysis-a hypothesis that requires further studies to verify.Overall, we postulate that the substrate and J-domain dependent allosteric pathways converge and that this convergence allows synchronization-only upon simultaneous interaction of substrate and J-domain is ATPase activity effectively stimulated, driving transition to the undocked state and substrate capture.
This coordinated mechanism also impacts the larger picture of Hsp70 allosteric dynamics.All Hsp70s studied share major characteristic transient conformational shifts-for example, undergoing repeated transitions involving undocking from the NBD of both SBDβ and SBDα.However, recent single molecule experiments indicate that the frequency of such transitions differs among members of the Hsp70 family [7-9, 18, 37-39].Interestingly, eukaryotic Hsp70s such as BiP of the endoplasmic reticulum [8,38] or Hsp70 and Hsc70 [39] of the cytosol, display conformational dynamics markedly different from DnaK and mitochondrial Hsp70.They have D SBD replaced by asparagine, pointing to the importance of the R NBD centered interactions for the allostery of Hsp70s.
Though no model for the mechanistic role of the D HPD -R NBD interaction per se has been previously put forth, the functional connection between these residues is well supported by genetic and biochemical data.Alteration of either R NBD or residues interacting with it (D SBD and D Lk ) [10,11,14,17,36,40] disrupt allosteric communication between domains, resulting in a somewhat increased basal activity and ineffective stimulation by its JDP partner, DnaJ in the case of DnaK [14] and Hsc20 in the case of Ssc1 (S4 Fig) .Substitution of R NBD (R167H) of DnaK suppresses nonfunctional D HPD (D35N) variant of DnaJ [20].However, we do not mean to say that this role of the D HPD is the only mechanistic effect of J-domain binding to Hsp70.For example, the invariant histidine and/or proline of the HPD may play specific roles, and, as previously suggested [41], conformational strain induced in the J-domain upon binding to Hsp70 might be important as well.However, there is no question that R NBD is key and could be considered as an allosteric perturbation site for the J-domain action.

Recognition step
At first glance, the sequence diversity of the interface of J-domain-Hsp70 complexes, but maintenance of overall positioning, and thus, function, seems surprising.It raises the question of how this diversity evolved without loss of function.Our analyses suggest that functional Jdomain/Hsp70 interactions have been maintained over time via coevolution of residues that form the binding interface.The nature of the Hsc20/Ssq1 interface is instructive in this regard.
Each key residue on one side of the large interface contacts more than one residue on the opposite side.Such a multifaceted interaction network is highly flexible.Over time residues can be added to or deleted from the interface without ever losing the ability to maintain functional interaction between the coevolving partners, enabling modulation of both the specificity and the strength of the interaction.Coevolution could also explain how a JDP might have been able to switch Hsp70 partnership during evolution, such as occurred in yeast mitochondriaswitching of Hsc20 from interacting with a general multifunctional Hsp70 (Ssc1) to a specialized Hsp70 (Ssq1) [42].A JDP that interacts with a given Hsp70 partner could gain the ability to interact with a "new" Hsp70 (i.e.product of gene duplication), upon emergence of mutation (s) that facilitate a new interaction, without affecting the initial one, at least during the period of transition [43].
What effect might this diversity have had on the functionality of Hsp70 chaperone systems more generally?We suggest that it has allowed evolution of complex Hsp70/JDP networks, both in situations in which a single type of Hsp70 interacts with many different JDPs and in which more than one Hsp70 coexist in a cellular compartment.In the first situation, replacement could well modulate the strength of a JDP/Hsp70 interaction, explaining why some JDPs have "easier access" to their partner Hsp70 than others.Here, electrostatic interactions could play a particularly important role.A JDP with a highly charged J-domain binding face would be able to recognize its Hsp70 partner at a longer distance, giving it priority over others in forming a productive interaction.In the second situation, the diversity serves to "insulate" Hsp70 networks, each Hsp70 interacting effectively only with its own subset of JDPs.For example, it could explain how yeast cytosolic Hsp70s Ssa and Ssb, which share common ancestry, evolved into independent networks, each with its own JDP partners [44,45].
In the two examples studied here, the predominant J-domain interactions are via helix II.However, their binding interfaces differ-the specialized Hsc20/Ssq1 interface involving mostly electrostatic interactions and the DnaJ JD /DnaK interface being a mosaic of electrostatic and hydrophobic interactions [22,41,46].However, across JDP/Hsp70s partnerships the interfaces may be even more varied.Interestingly, in two cases, Helix III of the J-domain was found to be very important for interaction with its Hsp70 partner-auxilin [35,47] and polyomavirus T Antigen [48].As only a few J-domain/Hsp70 interfaces have been studied, the breath of sequence variability, that supports proper positioning of the HPD sequence, the key to functional partnership, may be larger still.

Concluding remarks
Overall our data suggests a two-step process of J-domain function: initial binding that forms similar poses via highly variable residue-residue contacts, followed by a mechanistic step in which the invariant HPD perturbs the Hsp70 NBD-SBDβ,linker interface to foster a key allosteric transition.This separation could have future importance beyond a basic understanding of molecular chaperone function.It has long been appreciated that Hsp70 systems impinge on many biological functions and thus, not surprisingly, connected to many pathologic conditions raging from cancer, metabolic diseases to aging [49,50].As it has become appreciated that many Hsp70s partner with multiple JDPs [1,45], the idea of targeting a JDP, rather than Hsp70 itself, has emerged [51].Our findings showing a clear separation between recognition and mechanistic function provide a possible avenue to approach such intervention-modulation of the recognition step without affecting the mechanistic step, thus, manipulating the specificity or strength of particular J-domain/Hsp70 interactions.

Ssq1 homology model
Ssq1 model was built with I-TASSER [58] using the ATP-bound structure of DnaK (PDB ID 4B9Q) as a template with C-score 0.34.The Mg-ATP complex and structural water molecules were placed in the conserved ATP binding pocket of Ssq1 according to the DnaK structure (PDB ID 4B9Q).The obtained model of Ssq1 was refined by 3-μs-long MD simulation with position restraints applied to Mg-ATP and structural water during the first 1 μs to ensure a proper relaxation of the ATP-binding pocket.

Model of Hsc20-Ssq1 complex
The structures of Hsc20 (PDB ID 3UO3) and Ssq1 (homology model) used for molecular docking were first subjected to unbiased MD simulations to sample global protein structural dynamics.The Hsc20 simulation system was composed of a single protein in a 9.2 nm × 9.2 nm × 9.2 nm dodecahedron box filled with ~17000 water molecules, while the Ssq1 simulation system was composed of a single ATP-bound protein solvated with ~48,500 water molecules in a 13 nm × 13 nm × 13 nm dodecahedron box.Both systems were simulated using CHARMM36 [54] for 10 and 2 μs, respectively, and independently, using AMBER99SB-ILDN [59] for 5 and 2 μs, respectively, to test for force-field dependence.The obtained MD trajectories were separately clustered with g_cluster using a 0.35 nm RMSD cut off for Hsc20 and a 0.7 nm RMSD cut off for Ssq1 (to account for large conformational fluctuations of the sub-domain SBDα).Centroids of each of the 12 and 7 clusters identified for Hsc20 and Ssq1, respectively, were treated as representative structures of the proteins and paired together in 84 independent docking runs using ClusPro [60].As a result, 8,973 models were obtained.For each of the 84 docking runs, we selected the four best scoring models in each of the following criteria: a) binding energy estimation, b) minimal distance between the J-domain helix II and the NBD, c) minimal distance between the J-domain helix III and the NBD, d) minimal distance between the J-domain and the NBD.Redundancy among 336 selected models was removed with g_cluster, which resulted in 33 distinct models of the Hsc20-Ssq1 complex.Next, each of these models was placed in a 16 nm × 16 nm × 16 nm dodecahedron box and solvated with ~85000 water molecules.The resulting 33 systems were energyminimized and the protein side chains were relaxed with 50 ns-long MD simulations, during which the protein backbone atoms were restrained.Finally, the systems were simulated by unbiased MD for a minimum of 500 ns (total time for 33 simulated systems was 21 μs).The resulting trajectories were clustered based on heavy atoms of helices II and III of the Hsc20 J-domain after superimposing the Hsc20-Ssq1 complex structures using all Cα atoms of the conformationally rigid NBD domain.This was done using g_cluster with a single linkage method and a RMSD cutoff of 0.18 nm to obtain a fine characterization of the bound-state ensemble.

MD simulations of DnaK and the DnaJ JD -DnaK structures
The systems were composed of either DnaK (PDB ID 4B9Q) or DnaK in complex with the Jdomain of DnaJ (PDB ID 5NRO) placed in 13 nm × 13 nm × 13 nm dodecahedron box and solvated with approx.52,000 water molecules.Crystal water molecules coordinated to Mg 2+ ion in Mg-ATP were kept.Both systems were simulated for 1 μs and the final structures from these simulations were used as initial configurations for the subsequent free energy computations.

Energy landscape of intramolecular contacts with the conserved R NBD upon J-domain binding to Hsp70
The free energy landscapes describing the formation of intramolecular ion pairs by the critical R NBD (R207 in Ssq1 and R167 in DnaK) were computed with well-tempered multiple walkers 2D-metadynamics [61].In order to examine the effect of the J-domain binding on the R NBD interactions, the simulation systems contained either Hsc20-Ssq1 or DnaJ JD -DnaK complexes or individual proteins, Ssq1 or DnaK, respectively, placed in a 11.8 nm × 11.8 nm × 11.8 nm dodecahedron box filled with ~34,000 water molecules.To make the system significantly smaller and thus achieve better convergence of the free energy maps, the lid sub-domain, SBDα, was truncated by removing residues 570-657 of Ssq1 and 535-602 of DnaK.As two coordinates for 2D-metadynamics, we used the center of mass (COM) distances between the guanidine group of the R NBD residue and the carboxylic groups of the partner aspartate residues at the SBDβ (D517 in Ssq1 or D481 in DnaK) and at the interdomain linker (D429 in Ssq1 or D481 in DnaK).Biasing potential was added every 10 ps with the bias factor of 15, a height of 0.1 kJ and a width of 0.04 nm.The sampled range of both coordinates was restricted to a 0.3-1.9nm interval by applying one-sided harmonic potentials with a force constant of 3000 kJ nm -2 .For both systems, 10 independent 1-μs-long walker-simulations were carried out.The final free energy maps were averaged over 15 individual maps computed in 25-ns intervals from the last 350 ns of the simulations to average out free energy fluctuations [62].
The free energy profiles for the transition of the conserved arginine residue from intramolecular ion pairs with SBDβ (D517 for Ssq1, D481 for DnaK) and the linker (D429 for Ssq1, D481 for DnaK) to an intermolecular contact with the aspartate in the HPD (D50 for Hsc20, D35 for DnaJ) were computed using well-tempered multiple walkers metadynamics.The sampled coordinate was defined as the difference between the two distances: i) the COM distance between the arginine (R NBD ) guanidine group and the carboxylic groups of the aspartate residues in the linker (D Lk ) and in the SBDβ (D SBD ) and ii) the COM distance between the R NBD guanidine group and the carboxylic group of the D HPD (see Fig 2).To speed up convergence, both of these component distance coordinates were restricted to up to 1.5 nm by applying one-sided harmonic potentials with a force constant of 3000 kJ nm -2 .The bias potential was incremented every 10 ps by adding Gaussian-shaped potentials with height and width of 0.06 kJ and of 0.04 nm, respectively, and a bias factor of 8.For both systems, 10 walkers of 400 ns each were simulated at the same time.The final free energy profiles were averaged over 8 individual profiles computed in 20-ns intervals from the last 140 ns to limit free energy fluctuations.

Conformational transition of Hsp70 upon J-domain binding
The free energy profiles describing the compactness of the NBD, linker and SBDβ interface in Ssq1 and DnaK in the presence and absence of J-domain were computed with replicaexchange umbrella sampling (REUS).The simulations of the protein complexes (Hsc20-Ssq1 and DnaJ JD -DnaK) and free proteins (Ssq1 and DnaK) were performed in a 11.8 nm × 11.8 nm × 11.8 nm dodecahedron box filled with ~36,000 water molecules.The reaction coordinate was defined as the COM distance between the Cα atoms of SBDβ (residues 434-541 in Ssq1 and 395-505 in DnaK) and the Cα atoms of the NBD region at the interface with SBDβ (residues 112-143 and 186-262 in Ssq1, residues 75-104 and 146-226 in DnaK).To accelerate the free energy convergence, the lid sub-domain (SBDα) was truncated (residues 542-657 in Ssq1 and 506-602 in DnaK).For all systems, the initial configurations for REUS simulations were alternately selected from two independent 200 ns steered-MD simulations.In these simulations, the distance between SBDβ and NBD was gradually increased using a moving one-sided harmonic potential with a force constant of 2500 kJ nm -2 , applied to either the COM distance or the minimal distance between the Cα atoms of SBDβ and the NBD interfacial residues.The reaction coordinate in the 2.5-2.9 nm range was divided into 5 equally spaced windows in which the system was restrained using a harmonic potential with a spring constant of 2500 kJ nm -2 .The same approach was used to study the effect of the R NBD ->D substitutions (Ssq1 (R207D) and DnaK(R167D)) on the relative arrangement of NBD and SBDβ.To perform these simulations, the R NBD residue was replaced by aspartate in all initial structures.For each of the REUS windows at least 800 ns long simulation were performed.Exchanges between neighboring umbrella sampling windows were attempted every 2 ps.The first 100 ns of the trajectory was discarded and free energy profiles were computed using WHAM [63] with the Monte Carlo bootstrap method to estimate uncertainties of the free energy.The number of contacts between the domains (NBD/SBDβ or NBD/linker) is defined as the number of nonhydrogen atoms within a cutoff distance of 0.5 nm.To compute the number of contacts as a function of the COM distances between the domains, REUS trajectories were reweighted using a factor of exp((U i (r)−F i )/kT), where U i (r) denotes the value of the biasing potential for a given frame and F i denotes the free energy added to the i-th US window by the applied bias.

Hsc20-Ssq1 interaction analysis
All inter-protein interactions were analyzed using a 10.5 μs MD trajectory of the Hsc20-Ssq1 complex representing a dominant bound state.(i) Ion pairs between Hsc20 and Ssq1 were identified using a cut-off of 0.5 nm for the distance between the centers of mass of the interacting charged groups [64].(ii) Hydrogen bonds were identified using a standard geometric criterion: a 0.35 nm cut-off for donor-acceptor distance and a 30˚cut-off for hydrogen-donoracceptor angle.(iii) The hydrophobic contact surface area was estimated using the g_sas tool from the GROMACS package [52].The hydrophobic surface buried in the Hsc20/Ssq1 interface was computed as the difference between the surface of hydrophobic residues exposed to water in the dissociated and bound state.(iv) Relative contributions to the Hsc20-Ssq1 binding free energy of individual residues were calculated with g_mmpbsa [65] using a single trajectory approach, non-linear Poisson-Boltzmann solver and the CHARMM charges and atomic radii.

Spontaneous binding of Hsc20 to Ssq1 and binding free energy
To study spontaneous binding between Hsc20 and Ssq1, the simulation dodecahedron box of size 16.2 nm × 16.2 nm × 16.2 nm containing single copies of Hsc20 and ATP-bound Ssq1 solvated with approx.96500 water molecules was prepared.The proteins were initially placed such that the COM distance between positively charged residues of the Hsc20 J-domain helix II (R37, R38, R41) and negatively charged residues of Ssq1 NBD (D246, E248, D249, E253, D362, D364) was equal to 3 nm.For this system, five independent unbiased MD simulations were performed which were terminated after 60 ns as this time was sufficient for the native complex to spontaneously form in 2 out of 5 replicas.The free energy profiles for Ssq1 interaction with Hsc20 and with Hsc20(R37A,R41A) substitution variant were computed using replica-exchange umbrella sampling (REUS) [66].As a reaction coordinate for REUS simulations, we used the COM distance between the respective binding sites corresponding to residues 34-44 of Hsc20 (Helix II of the J-domain) and residues 244-255 and 425-429 of Ssq1 (NBD).The harmonic potential was used to restrain the system along the reaction coordinate (spring constants and centers of the potential are summarized in S4 Table ) and the initial configurations were taken from the spontaneous binding trajectories.The same initial configurations were used for the Hsc20(R37A,R41A) variant obtained by substituting R37 and R41 to alanine in all initial frames.For each of the REUS windows 400 ns long simulations were performed.Exchanges between neighboring windows were attempted every 2 ps and accepted by the Metropolis criterion.The free energy profiles were determined from the last 350 ns of thus obtained trajectories using the standard weighted histogram analysis method (WHAM) [63] and the Monte Carlo bootstrap method to estimate uncertainties of the free energy differences.
To infer Hsp70 phylogeny, a multiple sequence alignment was converted into a Hidden Markov Model [71] using hmmbuild program from the HMMER package.Forward-backward algorithm was used to compute a posterior probability (pp) for each site representing the degree of confidence in each position (residue or gap) of the alignment for each sequence.Amino acid positions with pp <0.7 were removed from the multiple sequence alignment [72].1,000 maximum likelihood (ML) searches were performed using RAxML v8.2.10 [73] with 1000 rapid bootstrap replicates, under the LG model of amino acid substitution and GAMMA model of rate heterogeneity with four discrete rate categories and the estimate of proportion of invariable sites (LG + I + G) [74], which was determined as the best-fit model by ProtTest v3.2 following Akaike Information Criterion [75].Hsp70s tree topology was used for all calculations.
Maximum likelihood implementation of the Coev model [32] was used to test coevolution between amino acid sequence positions homologous to those occupied by residues that were at a cutoff distance of 0.5 nm across the Hsc20/Ssq1 and DnaJ JD /DnaK binding interfaces.Coevolution analyses were based on the Hsp70 phylogeny, with Hsc20 orthologs paired with mitochondrial Hsp70 partners and DnaJ orthologs paired with DnaK partners.The fit of the Coev model in comparison to the 'independent' model was tested using Δ Akaike information criterion (ΔAIC = AIC independent -AIC Coev ).To establish the threshold for the ΔAIC value to be accepted as evidence for coevolution [32], we simulated sequence alignment based on the same phylogenetic tree as the original data but using the 'independent' model of sequence evolution.We used evolver software from the PAML 4 [76] for this simulation.The 99 th percentile of this expected ΔAIC distribution provided thresholds to consider the observed ΔAIC for Coev model to be accepted as coevolving with p<0.01 confidence level.
All other Hsc20 variants with a polyhistidine tag at the C-terminus were purified according to [24], except E. coli strain C41(DE3) was used for expression.Ssq1 variants were purified as described [77].Mge1 was purified as described [24].Isu1 from Chaetomium thermophilum (Isu1 Ct ) was purified as described [25] In all cases, protein concentrations, determined by using the Bradford (Bio-Rad) assay system with bovine serum albumin as a standard, are expressed as the concentration of monomers.

ATPase activity of Ssq1
The ATPase activity of Ssq1 variants defective in the Hsc20/Ssq1 interaction was measured as described by [24] with 0.5 μM Ssq1, 500 μM PVK-p peptide (LSLPPVKLHC) derived from the Isu1 sequence containing the PVK motif that interacts with substrate binding site of Ssq1 [78], 0.5 μM Mge1, and Hsc20 at the indicated concentrations in buffer A (40 mM Hepes-KOH, pH 7.5, 100 mM KCl, 1 mM dithiothreitol, 10 mM MgCl 2 , and 10% (v/v) glycerol).Reactions (15 μl) were initiated by the addition of ATP (2 μCi, DuPont NEG-003H, 3000 Ci/mmol) to a final concentration of 120 μM.Incubation was carried out at 25˚C, and the reaction was terminated after 15 min by the addition of 100 μl of 1M perchloric acid and 1mM sodium phosphate.After addition of 400 μl of 20 mM ammonium molybdate and 400 μl of isopropyl acetate, samples were mixed and the phases were separated by a short centrifugation.An aliquot of the organic phase (150 μl), containing the radioactive orthophosphate-molybdate complex, was removed and radioactivity was determined by liquid scintillation counting.Control reactions lacking protein were included in all experiments.Values were plotted in Graph-PadPrism.The ATPase activity corresponding to the maximal stimulation (Vm) of the wildtype Hsc20 was set to 1.

Site specific photo-crosslinking
Reaction mixtures (20 μl) containing 20 μM Hsc20 Bpa-containing variants and 5 μM Ssq1 variants were prepared in the reaction buffer (40 mM Hepes-KOH, pH 7.5, 100 mM KCl, 1 mM dithiothreitol, 10 mM MgCl 2 , and 10% (v/v) glycerol).Reactions were initiated by the addition of ATP to a final concentration of 2 mM and incubated at 25˚C for 15 min.Reactions were then exposed to 254 nm UV light (CL-1000 Ultraviolet Crosslinker) for 7 minutes.Control reactions were not exposed to UV light.Next, 5 μl of 4-fold concentrated Laemmli sample buffer (125 mM Tris-HCl, pH 6.8, 5% sodium dodecyl sulfate, 10% 2-mercaptoethanol, 20% (v/v) glycerol) was added to the reaction mixes.After 10 min incubation at 100˚C, the total reaction volume was loaded onto SDS-PAGE and visualized by Coomassie blue stain.Every experiment was replicated three times.

Identification of photo-crosslinking products using mass spectroscopy (MS)
Gel pieces were dried with acetonitrile and subjected to reduction with 10 mM dithiothreitol in 100mM NH 4 HCO 3 for 30 min at 57˚C.Cysteines were then alkylated with 0.5 M iodoacetamide in 100 mM NH 4 HCO 3 (45 minutes in a dark room at room temperature) and proteins were digested overnight with 10 ng/μl trypsin (CAT NO V5280, Promega) in 25 mM NH 4 HCO 3 [80].Digestion was stopped by adding trifluoroacetic acid at a final concentration of 0.1%.The mixture was centrifuged at 4˚C, 14 000 g for 30 min, to remove any remaining solid.MS analysis was performed by LC-MS in the Laboratory of Mass Spectrometry (Institute of Biochemistry and Biophysics, Polish Academy of Sciences, Warsaw) using a nanoAcquity UPLC system (Waters) coupled to an Orbitrap Elite (Thermo Fisher Scientific).The mass spectrometer was operated in the data-dependent MS2 mode, and data were acquired in the m/z range of 300-2000.Peptides were separated by a 180 min linear gradient of 95% solution A (0.1% formic acid in water) to 35% solution B (acetonitrile and 0.1% formic acid).The measurement of each sample was preceded by three washing runs to avoid cross-contamination.Data were searched with Mascot [81], with the following parameters: enzyme: Trypsin, parent ions mass tolerance: 20 ppm, fragment ion mass tolerance: 0.1, missed cleavages: 1, fixed modifications: Carbamidomethyl (C), variable modifications: Oxidation (M).Protein identification was validated using a target-decoy search strategy.

Fig 1 .
Fig 1. Binding mode of J-domain-Hsp70 complexes.(A) (top left) The dominant state of the Hsc20-Ssq1 complex obtained by protein-protein docking followed by MD refinement.Structure of S. cerevisiae Hsc20 (PDB ID 3UO3) was used along with a homology model of ATP-bound Ssq1, based on the X-ray structure of ATP-bound DnaK (PDB ID 4B9Q).The Ssq1 model had the expected architecture of an ATP-bound Hsp70: Nucleotide binding domain (NBD, grey) interacting with α and β subdomains of the substrate binding domain (SBD, brown) and the linker (yellow), in a groove between NBD subdomains Ia (light grey) and IIa (dark grey).To model the Hsc20-Ssq1 complex, conformational fluctuations of Hsc20 and ATP-bound Ssq1 were first characterized using MD simulations and clustering (S1 Fig).The obtained distinct conformers were then combinatorially docked.Based on the energetic and geometric criteria described in Methods, a representative set of 33 of the resulting complexes was selected for further refinement with MD simulations, and then clustered according to the positioning of the J-domain with respect to Ssq1 (S2 Fig).The model shown represents the dominant state with a population of 56%: J-domain of Hsc20 (cyan) interacts with Ssq1 at the NBD/SBDβ, linker interface.HPD sequence (red) of the J-domain is positioned near the conserved R207 (blue) of the NBD.(top middle) View of the Hsc20-Ssq1

Fig 2 .
Fig 2. D HPD of J-domain displaces interdomain contacts formed by R NBD with the D Lk and D SBD .(a) Free energy landscape governing the formation of Ssq1/DnaK R NBD interdomain contacts with D Lk and D SBD (top) Schematic depicting of the three Hsp70 residues of the R NBD interaction network; R NBD of subdomain Ia (blue dot) can interact with D LK of the linker (yellow dot) and D SBD of the SBDβ (brown dot).The four possible interaction states of R NBD are illustrated at right with the color of the rectangles serving as a key for the free energy landscapes at bottom.NBD subdomains Ia (light grey) and IIa (dark grey), SBD subdomains α and β (brown), linker (yellow).R NBD = R207/R167 in Ssq1/DnaK; D Lk = D429/D393 in Ssq1/DnaK; D SBD = D517/D481 in Ssq1/DnaK.(bottom) Free energy landscapes governing the formation of these interaction states for Ssq1 and DnaK alone (left panels) and for Hsc20-Ssq1 and DnaJ JD -DnaK complexes (right panels) with distances between RNBD-D Lk and between R NBD -D SBD used as relevant coordinates.Integrated populations of these interaction states are shown by pie charts (see S1 Table for data with standard errors).(b) Free energy profile of the transition of R NBD from intramolecular interactions with D Lk and D SBD residues to intermolecular contact with D HPD of the J-domain HPD loop when kept near NBD subdomain Ia.Distance difference Δd = d 1 -d 2 is used as a relevant coordinate.d 1 and d 2 are the distances between R NBD and D Lk /D SBD of Hsp70 (d1) or D HPD of the J-domain (d 2 ).Conformations in which R NBD interacts with D HPD at the expense of interaction with D Lk and/or D SBD (representative structures on the right, corresponding to the free energy minimum at 0.4-0.5 nm) are energetically more favorable than conformations with R NBD interacting simultaneously with D Lk and D SBD (representative structures on the left).
Fig 2. D HPD of J-domain displaces interdomain contacts formed by R NBD with the D Lk and D SBD .(a) Free energy landscape governing the formation of Ssq1/DnaK R NBD interdomain contacts with D Lk and D SBD (top) Schematic depicting of the three Hsp70 residues of the R NBD interaction network; R NBD of subdomain Ia (blue dot) can interact with D LK of the linker (yellow dot) and D SBD of the SBDβ (brown dot).The four possible interaction states of R NBD are illustrated at right with the color of the rectangles serving as a key for the free energy landscapes at bottom.NBD subdomains Ia (light grey) and IIa (dark grey), SBD subdomains α and β (brown), linker (yellow).R NBD = R207/R167 in Ssq1/DnaK; D Lk = D429/D393 in Ssq1/DnaK; D SBD = D517/D481 in Ssq1/DnaK.(bottom) Free energy landscapes governing the formation of these interaction states for Ssq1 and DnaK alone (left panels) and for Hsc20-Ssq1 and DnaJ JD -DnaK complexes (right panels) with distances between RNBD-D Lk and between R NBD -D SBD used as relevant coordinates.Integrated populations of these interaction states are shown by pie charts (see S1 Table for data with standard errors).(b) Free energy profile of the transition of R NBD from intramolecular interactions with D Lk and D SBD residues to intermolecular contact with D HPD of the J-domain HPD loop when kept near NBD subdomain Ia.Distance difference Δd = d 1 -d 2 is used as a relevant coordinate.d 1 and d 2 are the distances between R NBD and D Lk /D SBD of Hsp70 (d1) or D HPD of the J-domain (d 2 ).Conformations in which R NBD interacts with D HPD at the expense of interaction with D Lk and/or D SBD (representative structures on the right, corresponding to the free energy minimum at 0.4-0.5 nm) are energetically more favorable than conformations with R NBD interacting simultaneously with D Lk and D SBD (representative structures on the left).https://doi.org/10.1371/journal.pcbi.1007913.g002

Fig 4 .
Fig 4. Electrostatic interface of Hsc20 J-domain and Ssq1.(a) Network of key polar interactions at the binding interface between the J-domain of Hsc20 (cyan) and Ssq1 subdomain IIa (grey), linker (yellow) and SBDβ (brown).The most stable ion pairs and hydrogen bonds identified in MD simulations are indicated with dashed lines.Relative contributions of individual residues to the Hsc20-Ssq1 complex stability are indicated by color scale; blue for J-domain residues and red for Ssq1 residues.Note: K70 JD does not form obvious contacts but contributes strongly to the complex stability.(b) Site specific crosslinking of Hsc20 to Ssq1.Purified Hsc20 Q46Bpa (left) and M51Bpa (right) variants were incubated with purified Ssq1 wild-type (WT) or alanine substitution variants (D246A, E253A or E248A, D249A) in the presence of ATP.After UV irradiation (+), or as a control no irradiation (-), reaction mixtures were separated by SDS-PAGE.Migrations of size standard, in kDa, are indicated; Q:J � and Q: J �� indicate positions of the Ssq1-Hsc20 crosslinks, which were identified using mass-spectroscopy (S9 Fig); we do not understand the nature of the more slowly migrating crosslink band (Q:J �� ) present in the M51Bpa variant sample.Quantification of three independent Bpa crosslinking experiments; statistically significant differences (t-student test, p<0.001) are indicated ( ��� ).(c) Stimulation of Ssq1 ATPase activity by Hsc20 was measured in the presence of 0.5 μM of Ssq1, 500 μM of PVK-peptide substrate, 0.5 μM of nucleotide exchange factor Mge1 and indicated concentrations of Hsc20.(top) Hsc20 WT and alanine substitution variants.(bottom) Ssq1 WT and alanine substitution variants.Error bars indicate standard deviation for three independent measurements.https://doi.org/10.1371/journal.pcbi.1007913.g004

Fig 5 .
Fig 5. Sequence variability of J-domain/Hsp70 binding interface.(a) Pairwise sequence alignments of the J-domain/Hsp70 interfacial regions of the J-domains Hsc20 and DnaJ's (top) and the Hsp70s Ssq1 and DnaK (bottom).Residues of each protein that are in contact across the Hsc20/Ssq1 and DnaJ JD /DnaK binding interfaces are in green; HPD motifs in purple.(b) Sequence variability of positions that form the J-domain/Hsp70 binding interfaces (those positions shown in green in panel a) for orthologs of Hsc20/Ssq1 from 67 Ascomycota species (top) and DnaJ/DnaK from 117 bacteria species (bottom).Coevolving positions with p>0.01 statistical support (S2Table, S3 Table, S16 Fig and S17Fig) are connected by dashed lines.Sequence logos represent the amino acid frequency of each interfacial position in the orthologs examined, with positively charged (blue), negatively charged (red) and uncharged (grey) residues.Position numbering is that for S. cerevisiae (top) and E. coli (bottom) sequences.(c) Fraction of interfacial Hsc20/Ssq1 and DnaJ JD /DnaK residues that are positively and negatively charged.Whiskers indicate standard deviation.
Fig 5. Sequence variability of J-domain/Hsp70 binding interface.(a) Pairwise sequence alignments of the J-domain/Hsp70 interfacial regions of the J-domains Hsc20 and DnaJ's (top) and the Hsp70s Ssq1 and DnaK (bottom).Residues of each protein that are in contact across the Hsc20/Ssq1 and DnaJ JD /DnaK binding interfaces are in green; HPD motifs in purple.(b) Sequence variability of positions that form the J-domain/Hsp70 binding interfaces (those positions shown in green in panel a) for orthologs of Hsc20/Ssq1 from 67 Ascomycota species (top) and DnaJ/DnaK from 117 bacteria species (bottom).Coevolving positions with p>0.01 statistical support (S2Table, S3 Table, S16 Fig and S17Fig) are connected by dashed lines.Sequence logos represent the amino acid frequency of each interfacial position in the orthologs examined, with positively charged (blue), negatively charged (red) and uncharged (grey) residues.Position numbering is that for S. cerevisiae (top) and E. coli (bottom) sequences.(c) Fraction of interfacial Hsc20/Ssq1 and DnaJ JD /DnaK residues that are positively and negatively charged.Whiskers indicate standard deviation.https://doi.org/10.1371/journal.pcbi.1007913.g005 Fig. Representative structures of Ssq1 (a) and Hsc20 (b) from our unbiased molecular dynamics simulations carried out using the CHARMM36 (left) and AMBER99SB-ILDN (right) force fields.Ssq1 and Hsc20 were simulated for 2 and 10 μs, respectively, with CHARMM, and for 2 and 5 μs, respectively, with AMBER.The obtained trajectories were subject to cluster analysis with an RMSD cutoff of 0.7 nm (Ssq1) and 0.35 nm (Hsc20).The centroids of each cluster were superimposed and are shown in different colors.NBD, SBDa and SBDb domains (in Ssq1) and J-domain and Isu1 binding domain (in Hsc20) are indicated.(PDF) S2 Fig.(a) Five most populated binding poses of Hsc20 along with their percentage contributions to the bound-state ensemble of the Hsc20-Ssq1 complex, obtained by the combined docking/MD simulation approach.Helices II and III are shown in cyan and the rest of Hsc20 in transparent cyan.SBD is in light brown, Ia and IIa subdomains of NBD are in light and dark grey, respectively, linker in yellow.(b) Time evolution of the average RMSD of helices II and III of the J-domain with respect to the center of the most populated binding pose above (the dominant bound state).Averaging was done over all 18 independent MD trajectories that were found to visit the dominant bound state at least once over the course of 500 ns.The shaded area shows the standard deviation of the RMSD.(PDF) S3 Fig. Contact maps of Ssq1 R207 and DnaK R167 calculated from 2D metadynamics trajectories in absence (grey bars) and presence (cyan bars) of Hsc20 and DnaJ J-domain, respectively.The contact is defined with a cutoff distance of 0.5 nm between non-hydrogen atoms, interdomain linker (yellow); NBD subdomains Ia and IIa in light and dark grey, respectively.(c) Structural alignment of Hsc20-Ssq1 complex obtained by spontaneous binding MD simulations and the dominant bound state obtained by molecular docking/MD simulations.The heavy-atom RMSD for J-domains after aligning NBDs equals to 2.5 Å.(d) Free energy profiles for Hsc20-Ssq1 binding as a function of separation distance, d, between J-domain of Hsc20 and NBD subdomain IIa of Ssq1.Binding of Hsc20 wild-type (WT) (blue line), binding of Hsc20(R37A,R41A) variant (red line); shaded areas show the standard error.The vertical dashed line, demarcating the bund and unbound states, indicates the distance beyond which all specific interactions between the binding partners are lost.Representative structures are shown in circles: unbound state (left), bound state-corresponding to the energy well (right).(PDF) S12 Fig. Sequence variability of positions that form the Hsc20/Ssq1 binding interfaces (positions shown in green in Fig 5A) for orthologs of Hsc20 (left) and orthologs of Ssq1 (right) from the following taxonomic groups: Saccharomycotina (true yeasts), Pezizomycotina (p.yeasts), Taphiromycotina (t.yeasts), Basidiomycota (club fungi), Animals (animals) and Plants (plants).Sequence logos represent the amino acid frequency of each interfacial position in the orthologs examined, with positively charged (blue), negatively charged (red) and uncharged (grey) residues.Position numbering is that for S. cerevisiae.Phylogenetic relationships among taxonomic groups are depicted on the left.(PDF) S13 Fig. Sequence variability of positions that form the DnaJ/DnaK binding interfaces (positions shown in green in Fig 5A) for orthologs of DnaJ (top) and orthologs of DnaK (bottom) from bacteria and from mitochondria for the following taxonomic groups: Bacteria, Saccharomycotina (true yeasts), Pezizomycotina (p.yeasts), Taphiromycotina (t.yeasts), Basidiomycota (club fungi), Animals (animals) and Plants (plants).Sequence logos represent the amino acid frequency of each interfacial position in the orthologs examined, with positively charged (blue), negatively charged (red) and uncharged (grey) residues.Position numbering is that for E. coli.Phylogenetic relationships among taxonomic groups are depicted on the left.(PDF) S14 Fig. Maximum likelihood phylogeny of mitochondrial orthologs of Ssq1 from 67 Ascomycota species used in the sequence variability and co-evolution analyses.Major clades are indicated.(PDF) S15 Fig. Maximum likelihood phylogeny of 117 bacterial DnaK orthologs used in the sequence variability and coevolution analyses.Major clades are indicated.(PDF) S16 Fig. ΔAIC (AIC independent -AIC Coev ) distribution calculated for each pair of positions derived from amino acid sequence of length 33 simulated under the independent evolution model using mitochondrial Hsp70 phylogeny (S12 Fig).The co-evolution acceptance threshold with p>0.01 confidence was established as ΔAIC = 7.38 (red dashed line).(PDF) S17 Fig. ΔAIC (AIC independent -AIC Coev ) distribution calculated for each pair of positions derived from amino acid sequence of length 54 simulated under the independent evolution model using DnaK phylogeny (S13 Fig).The co-evolution acceptance threshold with p>0.01 confidence was established as ΔAIC = 9.473 (red dashed line).(PDF)