Sequence-specific Long Range Networks in PSD-95/Discs Large/ZO-1 (PDZ) Domains Tune Their Binding Selectivity*

Protein-protein interactions mediated by modular protein domains are critical for cell scaffolding, differentiation, signaling, and ultimately, evolution. Given the vast number of ligands competing for binding to a limited number of domain families, it is often puzzling how specificity can be achieved. Selectivity may be modulated by intradomain allostery, whereby a remote residue is energetically connected to the functional binding site via side chain or backbone interactions. Whereas several energetic pathways, which could mediate intradomain allostery, have been predicted in modular protein domains, there is a paucity of experimental data to validate their existence and roles. Here, we have identified such functional energetic networks in one of the most common protein-protein interaction modules, the PDZ domain. We used double mutant cycles involving site-directed mutagenesis of both the PDZ domain and the peptide ligand, in conjunction with kinetics to capture the fine energetic details of the networks involved in peptide recognition. We performed the analysis on two homologous PDZ-ligand complexes and found that the energetically coupled residues differ for these two complexes. This result demonstrates that amino acid sequence rather than topology dictates the allosteric pathways. Furthermore, our data support a mechanism whereby the whole domain and not only the binding pocket is optimized for a specific ligand. Such cross-talk between binding sites and remote residues may be used to fine tune target selectivity.

Several studies on protein-protein interaction modules (1), such as SH2, SH3, and PDZ 5 domains, have highlighted the apparent problem of how selectivity is achieved in a cellular context (2)(3)(4). For example, the synapse contains several PDZ proteins (5) with apparently overlapping specificity, and yet they have distinct roles in scaffolding and signaling (6). One mechanism that could modulate specificity is intradomain allostery (7). Multisubunit proteins such as hemoglobin are well known for their allosteric behavior whereby the oxygen binding is cooperative because ligation of one subunit will modulate the affinity in the other subunits of the tetramer (8). But allostery is also present within monomeric proteins (9), and it may be invoked even in the absence of a detectable conformational change (10 -15).
Previous studies on monomeric globular proteins have predicted allosteric pathways on an amino acid residue level based on analysis of co-evolution (16,17), molecular dynamics simulations (18), or NMR-constrained molecular dynamics (19). Although the PDZ domain family has served as a prime model system for these studies, experimental data that detect intramolecular allosteric pathways are scarce and sometimes conflicting (16,20,21).
In protein folding studies, it proved illuminating to compare homologous proteins to unravel basic determinants of the underlying mechanism (22). Here, we employ the same strategy to assess the structural pattern of intradomain allostery by looking at the coupling free energy (⌬⌬⌬G C ) upon binding of peptide ligands to two different members of the PDZ domain family. We used a double mutant cycle analysis (23) to quantify energetic coupling between different residues in the binding reaction (see Fig. 1).
The protein domains chosen for the study, protein tyrosine phosphatase Basophil-like PDZ2 (hereafter called PDZ2) and PSD-95 PDZ3 (PDZ3), are perfectly suited for this type of analysis on allosteric networks because they have been extensively investigated by experiments directed to folding and binding (21, 24 -26). Importantly, computational studies on these PDZ domains predicted the presence of distinct allosteric networks (16,18,19,27,28). In this study, we present experimental evidence for such networks.

Mutagenesis, Expression, and Purification of Recombinant
PDZ2 and PDZ3 Domains-Site-directed mutants were prepared by inverted PCR using the cDNAs encoding pseudo-wild type PDZ2 and pseudo-wild type PDZ3, respectively, as templates. These cDNAs code for a Trp residue in homologous positions for the two domains, namely Y43W and F337W, respectively, as probe for binding (20,24,25,29). PDZ2 and PDZ3 were expressed in Escherichia coli and purified as described previously (30,31). Protein concentrations were determined by absorbance measurements using calculated extinction coefficients from the amino acid composition (PDZ3, ⑀ 280 ϭ 8500 M Ϫ1 cm Ϫ1 ; PDZ2, ⑀ 280 ϭ 5500 M Ϫ1 cm Ϫ1 ). The numbering of residues in this article is in accordance with NMR and crystal structures of PDZ2 and PDZ3: Protein Data Bank codes 1VJ6 (25) and 1BE9 (32), respectively.
Ligand Binding Kinetics-Single-mixing kinetic binding experiments were carried out on SX-18MV and SX-20MV stopped-flow instruments (Applied Photophysics, Leatherhead, UK) by following the change in fluorescence emission upon binding. Excitation was at 280 nm, and emission was measured by using either a 330 nm band-pass filter or a 435 nm cut-off filter. We thus used fluorescence resonance energy transfer between the tryptophan in the PDZ domain and a dansyl group covalently bound to the N terminus of the specific peptide. We have shown previously that this dansyl does not affect the binding kinetics of hexapeptides (20,31,33), and any small effects would cancel out in the calculation of ⌬⌬G for wild type and modified peptides.
All mutants were subjected to kinetic binding measurements using three different peptide ligands: i) the wild type peptides, EQVSAV for PDZ2 as described in Gianni et al. (25) and YKQTSV derived from the C terminus of the protein CRIPT (34) for PDZ3; ii) the wild type peptides containing a Val 0 3 Abu mutation for both PDZ domains and iii) a Ser Ϫ2 3 Thr or Thr Ϫ2 3 Ser mutation in the peptides for PDZ2 and PDZ3, respectively. Abu is an ␣-aminobutyric acid, i.e. a Val with one ␥-methyl group deleted. Peptides were purchased from GL Biochem (Shanghai, China) and JPT Peptide Technologies (Berlin, Germany). The concentration of peptide was determined by measuring the dansyl absorbance at 330 nm using ⑀ 330 ϭ 4,300 M Ϫ1 cm Ϫ1 .
All measurements on PDZ2 were performed using 2 M protein concentration in 50 mM phosphate buffer (pH 7.2) and 0.4 M Na 2 SO 4 at 4°C. Experiments with PDZ3 were performed using 3 M protein concentration in 50 mM phosphate (pH 7.45) at 10°C. In all cases, the time course for binding was satisfactorily fitted by a single exponential function. The observed reaction was thus consistent with the simplest kinetic mechanism involving a bimolecular association rate constant and a monomolecular dissociation rate constant. Kinetic traces (supplemental Figs. S1A and S2A) were fitted to a single exponential equation using the software provided by Applied Photophysics to obtain the observed rate constant, k obs (Equation 1), where A is the fluorescence signal at time t, ⌬A EQ is the amplitude of the kinetic trace, and C is the fluorescence at the beginning of the reaction.
The k obs values were then plotted versus the peptide concentration and analyzed by either Equation 2, which is valid under pseudo-first-order conditions (PDZ2) or the general equation for a second-order bimolecular association (35) for PDZ3 (Equation 3).
In these equations, k on is the association (on-rate) constant, k off is the dissociation (off-rate) constant, and [A] 0 and n are the initial concentrations of PDZ and peptide, respectively.
Displacement Kinetics-The off-rate constant k off for PDZ3 was in many cases Ͻ5 s Ϫ1 and could therefore not be estimated with good enough accuracy from Equation 3. All k off values for PDZ3 were therefore determined in displacement reactions (33). PDZ3 and dansylated peptide were first preincubated (0.5-3 M PDZ and 0.25-3 M dansylated peptide, depending on the respective affinity between the particular pair of PDZ and peptide). The complex was mixed rapidly with an excess of unlabeled peptide. At sufficiently high concentration of unlabeled peptide the dissociation of dansylated peptide from the protein is the rate-limiting step (k off ) (supplemental Fig. S2, C and D). Kinetic traces were single exponential and fitted to Equation 1 to obtain observed rate constants (supplemental Fig.  S2C). These k obs values were plotted versus the concentration of the unlabeled peptide used to displace the complex (supplemental Fig. S2D). The data were analyzed by Equation 4 to estimate a value for k off . This equation is valid under pseudofirst-order conditions, i.e. at high concentrations of unlabeled peptide. But, it also serves as a good model-free equation at lower concentrations of unlabeled peptide, where K D is the equilibrium dissociation constant for the PDZ and the unlabeled peptide and k on Ј the apparent (first-order) on-rate constant for the labeled peptide and the PDZ domain at the chosen concentration of labeled peptide. Note that k off can be accurately estimated from this experiment as the asymptotic value of k obs at high concentration of unlabeled peptide. All data were fitted to Equations 2-4 using KaleidaGraph software (Synergy Software).
Calculation of Coupling Energies-Coupling energies were calculated as described (23). In brief, for a hypothetical wild type protein P-AB, the changes in free energy for single mutants at positions A and B can be calculated assuming the following (Equations 5-7).
P-A is the variant where B is deleted, P-B is the variant where A is deleted, and P is the variant where both A and B are deleted. By applying a thermodynamic cycle (Equation 5) to the four species P-AB, P-A, P-B, and P, it may be noted that if A and B do not interact energetically, the change in energy upon deleting A in P-AB (⌬⌬G P-AB3 P-B ) should equal the change in energy upon deleting A in P-A (⌬⌬G P-A3 P ). Thus, a coupling energy ⌬⌬⌬G C may be calculated as the difference in ⌬⌬G for deleting A in P-AB and for deleting A in P-A.
A detectable ⌬⌬⌬G C 0 implies that the probed positions A and B interact energetically, and its value is a quantitative measure of the strength of this interaction.
Fersht and co-workers (23,36,37) initially developed the double mutant cycle methodology for probing intramolecular interactions. But, as shown by Schreiber and Fersht (38,39), it works equally well for intermolecular interactions, where the energetic interaction between pairs of residues on different molecules can be monitored directly. This first and most detailed analysis of a protein-protein interaction by double mutant cycles was on that between barnase and barstar.
In the case of intermolecular interactions, binding rate constants are measured for the four different species PA (wild type protein), P (single mutant protein where A is deleted), LB (wild type ligand), and L (single mutant ligand where B is deleted). The experimental observables are the association and dissociation rate constants for PA with LB, i.e. k on(PA-LB) and k off(PA-LB) , and to L, k on(PA-L) and k off(PA-L) , and the association and dissociation rate constants for the mutant protein P with LB, i.e. k on(P-LB) and k off(P-LB) , and to L, k on(P-L) and k off(P-L) .
By applying the double mutant cycle method (see Fig. 1), a coupling energy ⌬⌬⌬G C may be calculated as follows.
Two examples of calculation of ⌬⌬⌬G C are given in supplemental Figs. S1 and S2 for PDZ2 and PDZ3, respectively.

RESULTS
Allosteric systems are characterized by an energetic coupling between remote sites in a protein. Indeed, if allostery is at play, the existence of functional energetic networks between residues located in the ligand binding pocket or even the ligand itself, and those in other parts of the protein should be detectable. Pinpointing such networks is the focus of this study.
Double Mutant Cycles-A powerful approach to address quantitatively the energetic connectivity between remote residues in a protein is the double mutant cycle analysis (Fig. 1). This methodology has been used to capture the strength and the location of both intra-and intermolecular interactions in proteins (23,38).
The basic experimental parameters to be determined are the rate constants for the association and dissociation of a specific ligand to the known binding site on the protein. In a double mutant cycle analysis the kinetics of the wild type protein are compared with those of two single mutants at positions A and B and of the corresponding double mutant (both A and B being mutated). In our study, residue A belongs to the protein, the PDZ domain, and residue B to its specific peptide ligand.
We have studied 31 different mutants for PDZ2 and 31 mutants for PDZ3. The binding kinetics for each one of these mutants have been characterized for three different peptide ligands. The change in binding free energy of the wild type peptide LB, upon mutation of residue A, is quantified by ⌬⌬G PA-LB3 P-LB (Equation 9). Likewise, the change in binding free energy of the mutant peptide L, upon mutation of residue A, is quantified by ⌬⌬G PA-L3 P-L (Equation 10). If the two residues A and B are functionally independent in the binding reaction, the effect of mutating residue A in the protein would be identical in the binding reactions with wild type peptide LB and mutant peptide L, respectively. On the other hand, if the residues A and B are energetically coupled, ⌬⌬G PA-LB3 P-LB and ⌬⌬G PA-L3 P-L will not be identical and therefore a detectable coupling energy ⌬⌬⌬G C , according to Equation 11, will be observed.
Binding Kinetics of PDZ Mutants with Their Peptide Ligands-Side chains of core residues are likely candidates for transmitting allosteric signals through a protein domain (40). The two PDZ domains were extensively mutated such that all regions of the hydrophobic core and a few surface patches were probed by conservative mutations. All 62 mutants, although typically destabilized relative to the wild type, were all folded under the experimental conditions (24,29).
The binding rate constants were determined by fluorescence spectroscopy using the stopped-flow technique (supplemental Figs. S1 and S2). Association (k on ) and dissociation (k off ) rate constants with three different hexapeptide ligands (Fig. 2) were determined for each PDZ mutant (supplemental Table S1). These peptide ligands were as follows: (i) the wild type peptide corresponding to a previously identified natural ligand specific for each PDZ domain and (ii) peptides containing a Val 0 3 Abu mutation (deletion of a ␥-methyl group from Val 0 ) or (iii) a Ser Ϫ2 3 Thr or Thr Ϫ2 3 Ser mutation for PDZ2 and PDZ3, respectively (see Fig. 2 for numbering of the peptide).
We chose these peptide mutations because they target residues that interact with the binding pocket of PDZ domains. Furthermore, deletion of a methyl group is a relatively conservative mutation (41) that perturbs the affinity between PDZ and peptide to a perfect degree for analysis (usually 0.5-1 kcal mol Ϫ1 ) and for our stopped flow measurements (k off values between 0.1-300 s Ϫ1 ).
The wild type peptide for PDZ2 contains at the Ϫ2 position a Ser (Fig. 2) that could not be mutated to Ala because the ensuing very low affinity between PDZ and the ligand would jeopardize analysis (20); thus, the Ser Ϫ2 3 Thr mutation was chosen to probe this position. Importantly, PDZ2 is expected to display a poor selectivity for Ser versus Thr at position Ϫ2 of the ligand (42). Thus, in this case, the dissociation equilibrium constant K D for the PDZ2-peptide complex is essentially unchanged (supplemental Table S1).
Coupling Free Energies between PDZ and Peptide Side Chains-The coupling free energies between each mutated residue in the protein and the side chain of the two peptide side chains were calculated using the determined rate constants k on and k off (Table 1). We found that both PDZ domains display long range communication between the peptide and residues in the protein that are spatially distant from the binding site (up to 15 Å C ␣ -C ␣ distance), thereby supporting an allosteric mechanism.
The structural distribution of the coupling energies for the different residues can be seen in Figs. 2 and 3 for each PDZ domain. By comparing the two, it is immediately apparent that the wiring of the long range networks of coupling energies are different for the two PDZ domains, indicating distinct allosteric pathways. The overall coupling pattern for PDZ3 suggests a segment through the core with higher coupling free energies. This segment extends through the core at an angle compared with the two ␤-sheets of the PDZ-peptide complex. Furthermore, it appears to include the Thr Ϫ2 of the peptide (Fig. 2, B  and D). For PDZ2, on the other hand, the observed segment of higher coupling free energies goes from the loop between ␤2 and ␤3 (bottom of Fig. 2, A and C) and along the bound peptide toward Val 0 .
This interesting feature of specific coupling patterns is further demonstrated by the lack of correlation when the coupling   (Fig. 4A). This is in stark contrast to folding ⌽ values, which correlate strongly for these two protein domains, suggesting a common folding mechanism (24). It thus appears that whereas folding pathways are conserved within the PDZ domain family, their functional energetic networks are more variable. In other words, although native topology sculpts the free energy landscape for folding in a rather common way, the amino acid composition plays a key role in allostery. The sequence dependence of the allosteric pathways is illustrated by some of the positions with the highest coupling energy. For example, Val 29 in PDZ2 couples strongly to Val 0 in the peptide, whereas the corresponding Ile 327 in PDZ3 has a ⌬⌬⌬G C Ͻ0.2 kcal mol Ϫ1 . However, the same Ile 327 in PDZ3 has one of the highest ⌬⌬⌬G C for Thr Ϫ2 (0.48 kcal mol Ϫ1 ), whereas ⌬⌬⌬G C between the corresponding positions in PDZ2 is close to zero.
Nevertheless, a few positions display high coupling for both PDZ domains, for example position Leu 25 in PDZ2 and Leu 323 in PDZ3. This position was predicted as an allosteric network residue (18,19,28), and indeed, it couples with the C-terminal Val of the peptide (Table 1 and supplemental Fig. S1).
Similarly, a subset of residues couples to both the Val 0 and Ser/Thr Ϫ2 positions (e.g. Leu 25 , Thr 35 , Val 65 , Leu 73 , Leu 85 , and  Leu 94 in PDZ2 and Leu 323 , Phe 340 , Ile 341 , and Lys 380 in PDZ3; Table 1). Yet, no general correlation was observed between residues that couple to position Val 0 of the peptide and residues that couple to position Ser/Thr Ϫ2 (Fig. 4B). Interestingly, for PDZ2, Leu 25 , Val 65 , and Leu 94 have positive coupling energies with position Val 0 but negative coupling energies with Ser Ϫ2 .

DISCUSSION
The crowded and heterogeneous cellular environment demands the interactions between its numerous components to be highly specific, to avoid potentially harmful side reactions.
In fact, it is somewhat surprising that the cell employs effectively a limited number of distinct protein-protein interaction domains (2), such as the PDZ family, which is involved in the control of quite different cellular events. Despite the fact that this protein family displays a simple topology and a rather conserved binding pocket, the distinct roles of the various PDZ domains in cellular signaling require selectivity. Until recently, PDZ domains were grouped into three classes, depending on their consensus target sequence. A detailed analysis of the PDZ domains in the mouse proteome, however, suggested that separation into discrete classes might not be justified; instead, PDZ domains are evenly distributed throughout selectivity space, suggesting they have been optimized to increase the repertoire of specificities but maintaining minimal cross-reactivity (4). In this perspective, it is of interest to discuss the general significance of the coupling free energies that we obtained for the different peptides, which seem to imply a mechanism of optimization involving the whole domain rather than the binding site only.
The experimental data reported above show that all the mutations of PDZ2 that we have explored lead to a positive coupling energy for the Val 0 to Abu substitution (Fig. 3), although the absolute values are quite variable. That is, mutation of any core residue in the protein will affect the binding of the peptide such that the effect of deleting a methyl group from the native C-terminal Val 0 of the peptide will be smaller in the mutant as compared with the wild type PDZ2. This intriguing finding suggests that the entire PDZ scaffold has been under selective pressure to optimize the binding of the ligand by intradomain allosteric coupling. Sequence variation in the protein domain is therefore not neutral to peptide binding, strongly indicating that selectivity by the domain is not solely determined by the subset of residues directly involved in ligand binding. Additional support for this conclusion comes from the effect of the Ser 3 Thr mutation in the peptide for the binding with PDZ2. In this case, the rather even distribution of negative coupling energies (Figs. 3 and Fig. 5C) suggests that this position is not optimized for Ser relative to Thr (42), as also indicated by the nearly identical K D for wild type PDZ2 and the two peptide variants (supplemental Table S1).
Likewise, PDZ3 also displays mainly positive coupling energies. However, in this case, deletion of the methyl group of the Thr is more detrimental to the binding than the Val 3 Abu substitution (supplemental Table S1). Indeed, the Thr Ϫ2 3 Ser mutation displayed an even higher number of positive coupling energies than Val 0 3 Abu (Figs. 3 and 5). In conclusion, the high number of positive ⌬⌬⌬G C values suggests that PDZ3 is optimized for both Thr at position Ϫ2 and Val at position 0.
The coupling free energy between pairs of residues was found previously to depend strongly on the distance between the two residues (38,43). A similar conclusion was made for PDZ3 when coupling to the evolutionarily conserved residue His 372 was assessed, although with a limited data set (20). This behavior is not surprising because neighboring residues are expected to interact more effectively than distant ones. However, in the case of both PDZ2 and PDZ3, we did not observe a clear trend in the dependence of the coupling energies on distance, for neither of the peptide residues Val 0 and Ser/Thr Ϫ2 (Fig. 5). This result corroborates the hypothesis that sequence rather than topology governs the observed inter-residue communication. It is also in line with the observed distinct segments of higher coupling free energies for PDZ2 and PDZ3, respectively (Fig. 2).
In summary, we identified by double mutant cycles and binding kinetics, energetically coupled residues in two different homologous protein-ligand complexes. Our experimental results confirm the existence of long range networks suggested to exist in different proteins (17) and in particular predicted in the PDZ domain family (e.g. Refs. 16,18,19,27,28). Several predictions capture overall features and even details (18,19) of our present results. Our data may thus be used to benchmark the extensive available and forthcoming theoretical work on the PDZ system.
We suggest that distinct networks will be found in most proteins, but the extent of functional selectivity by a particular allosteric network will have to be assessed in each case. Nevertheless, these coupling networks exist and may be of advantage for a specific function and fixed by natural selection. Overall, our data suggest an attractive scenario whereby a subtle allosteric mechanism reduces cross-reactivity by optimizing ligand recognition using the whole domain. The mechanism involves distinct pathways that modulate specificities in different PDZ domains.