Noah Giacon,Ettore Lo Cascio,Darcy S. Davidson,Marcelo D. Polêto,Justin A. Lemkul,Valeria Pennacchietti,Livia Pagano,Carlotta Zamparelli,Angelo Toto,Alessandro Arcovito
Monomeric and dimeric states of human ZO1-PDZ2 are functional partners of the SARS-CoV-2 E protein
SARS-CoV-2 E protein,Human ZO1-PDZ2 domain,DRUDE force field,Molecular dynamics,Surface Plasmon Resonance,Folding upon binding
The Envelope (E) protein of SARS-CoV-2 plays a key role in virus maturation, assembly, and virulence mechanisms. The E protein is characterized by the presence of a PDZ-binding motif (PBM) at its C-terminus that allows it to interact with several PDZ-containing proteins in the intracellular environment. One of the main binding partners of the SARS-CoV-2 E protein is the PDZ2 domain of ZO1, a protein with a crucial role in the formation of epithelial and endothelial tight junctions (TJs). In this work, through a combination of analytical ultracentrifugation analysis and equilibrium and kinetic folding experiments, we show that ZO1-PDZ2 domain is able to fold in a monomeric state, an alternative form to the dimeric conformation that is reported to be functional in the cell for TJs assembly. Importantly, surface plasmon resonance (SPR) data indicate that the PDZ2 monomer is fully functional and capable of binding the C-terminal portion of the E protein of SARS-CoV-2, with a measured affinity in the micromolar range. Moreover, we present a detailed computational analysis of the complex between the C-terminal portion of E protein with ZO1-PDZ2, both in its monomeric conformation (computed as a high confidence AlphaFold2 model) and dimeric conformation (obtained from the Protein Data Bank), by using both polarizable and nonpolarizable simulations. Together, our results indicate both the monomeric and dimeric states of PDZ2 to be functional partners of the E protein, with similar binding mechanisms, and provide mechanistic and structural information about a fundamental interaction required for the replication of SARS-CoV-2.
Monomeric and dimeric states of human ZO1-PDZ2 are functional partners of the SARS-CoV-2 E protein

The Envelope (E) protein of SARS-CoV-2 plays a key role in virus maturation, assembly, and virulence mechanisms. The E protein is characterized by the presence of a PDZ-binding motif (PBM) at its C-terminus that allows it to interact with several PDZ-containing proteins in the intracellular environment. One of the main binding partners of the SARS-CoV-2 E protein is the PDZ2 domain of ZO1, a protein with a crucial role in the formation of epithelial and endothelial tight junctions (TJs). In this work, through a combination of analytical ultracentrifugation analysis and equilibrium and kinetic folding experiments, we show that ZO1-PDZ2 domain is able to fold in a monomeric state, an alternative form to the dimeric conformation that is reported to be functional in the cell for TJs assembly. Importantly, surface plasmon resonance (SPR) data indicate that the PDZ2 monomer is fully functional and capable of binding the C-terminal portion of the E protein of SARS-CoV-2, with a measured affinity in the micromolar range. Moreover, we present a detailed computational analysis of the complex between the C-terminal portion of E protein with ZO1-PDZ2, both in its monomeric conformation (computed as a high confidence AlphaFold2 model) and dimeric conformation (obtained from the Protein Data Bank), by using both polarizable and nonpolarizable simulations. Together, our results indicate both the monomeric and dimeric states of PDZ2 to be functional partners of the E protein, with similar binding mechanisms, and provide mechanistic and structural information about a fundamental interaction required for the replication of SARS-CoV-2.

Protein-protein interactions (PPIs) govern and regulate the vast majority of physiological and molecular pathways in the cell, representing key events for every aspect of cellular life. PPIs are usually mediated by globular domains that can bind and recognize so called short linear motifs (SLiMs) characterized by specific consensus sequences. Thus, it is not surprising that many pathogens, from viruses to bacteria, evolved to mimic those SLiMs in their proteins with the ultimate goal to acquire the ability to interact with host proteins, disrupt normal cell physiology and eventually lead to diseases. The Envelope (E) protein is one of the structural proteins encoded by the genome of SARS-CoV-2, the causative virus of Covid-19 pathology [1], [2]. Together with the Nucleocapsid protein (N), the Spike protein (S) and the Membrane protein (M), the E protein contributes to form a complete and functional viral structure [2], [3]. It is the shortest structural protein encoded by coronaviruses and consists of an integral membrane protein of 75 residues arranged in three domains: a short hydrophilic N-terminal domain (NTD), a large, α-helical hydrophobic transmembrane domain (TMD), and the largest hydrophilic C-terminal domain (CTD) [4]. Curiously, while the coronaviral envelope consists predominantly of M with only a small portion of E that is incorporated into the viral envelope of virions, the E protein is largely expressed in the host cell during virus replication, demonstrating an active role in virion trafficking, assembling and budding [5]. One of the structural key elements of the E protein that allows this complex set of functions is a PDZ-binding motif (PBM) contained in the CTD, that allows the viral protein E to interact with host cell proteins through one or more PDZ domains. PDZ domains (taking their name from PSD-95/Dlg/ZO1 proteins) are the most abundant protein-protein interaction modules of the human proteome, with more than 400 human proteins contain at least one PDZ domain in their sequence. PBMs are short, linear motifs of 4 residues usually located at the extreme C-terminus of PDZ-interacting proteins, although in some cases they can be found at intermediate positions along the polypeptide sequence. They can be classified based on their sequence, which determines and regulates binding specificity [6]. Notably, the binding capability of the C-terminal portion of the SARS-CoV-2 E protein with different PDZ-containing proteins involved in the formation of epithelial tight junctions (TJs) and in the regulation of cell-polarity processes has been reported [7], [8], [9], [10]. The interactions between the E protein and PDZ-containing proteins have been suggested to be the basis of the disruption of lung epithelial tissues occurring upon infection by coronaviruses [11], [12]. ZO1 is a PDZ-containing protein with a fundamental role in the formation of epithelial and endothelial TJs, multiprotein complexes that primarily serve to prevent the diffusion of ions and macromolecules through the paracellular environment between epithelial cells [13] and contribute to maintain the integrity of epithelial and endothelial tissues. ZO1 is a component of the cytosolic protein complex that interacts with transmembrane junctional proteins of TJs. The N-terminal portion of ZO1 is composed of three PDZ domains (PDZ1, PDZ2 and PDZ3), a SH3 domain and a Guanylate-Kinase Homology domain (GUK), while the C-terminal portion is variable due to alternative splicing in different tissues[14], [15]. The PDZ domains and the SH3 domain are the mediators of protein-protein interactions between ZO1 and different partners, regulating the spatial arrangement of the cytosolic protein complex of TJs [16], [17], [18], [19]). Interestingly, structural characterizations of ZO1 suggest the formation of homodimer complexes occurring through domain swapping at the level of PDZ2 domain [20]. In particular, the constitution of ZO1 homodimer has been proposed at the basis of its function, underlying the scaffolding and physiological properties of ZO1 [20]. In this work, through a combination of experimental and theoretical approaches, we characterized the binding between a peptide mimicking the SLiM in the C-terminal portion of the E protein from SARS-CoV-2 and the PDZ2 domain of ZO1 from molecular and structural perspectives. Specifically, we directly measured the binding affinity of PDZ2 for the E protein by surface plasmon resonance (SPR), finding an equilibrium dissociation constant in the order of micromolar range, remarkably lower than the nanomolar affinity that was previously measured for this complex [21]. Interestingly, we established PDZ2 to be monomeric in our experimental conditions, highlighting its ability to bind and recognize substrates as a monomer. By employing molecular dynamics (MD) simulations, we analyzed the structural determinants of the binding process, providing a molecular rationalization of the complex formation. Our data is discussed in light of available literature about the role of the SARS-CoV-2 E protein in Covid-19 pathogenicity and the underlying implications of its ability to bind PDZ2 domain of ZO1.

The construct encoding the PDZ2 domain of ZO1 protein was subcloned in a pET28b+ plasmid vector and then transformed in Escherichia coli BL21 (DE3) cells. Bacterial cells were grown in LB medium, containing 30 μg/mL of kanamycin, at 37 °C until OD600 = 0.7 − 0.8, and then protein expression was induced with 0.5 mM IPTG. After induction, cells were grown at 25 °C overnight and then collected by centrifugation. The pellet was resuspended in buffer containing 50 mM TrisHCl, 300 mM NaCl, 10 mM imidazole, pH 8.0, adding one antiprotease tablet (Complete EDTA-free, Roche), then sonicated and centrifuged. The soluble fraction from the bacterial lysate was loaded onto a Ni-charged HisTrap Chelating HP (GE Healthcare) column. The equilibration buffer was 50 mM TrisHCl, 300 mM NaCl, 10 mM imidazole, pH 8.0. Elution was performed with an ÄKTA-prime system, with a gradient of imidazole from 0 to 1 M. Fractions containing the protein were determined through SDS-PAGE. The buffer was exchanged to 10 mM Hepes, 150 mM NaCl, pH 7.4 by using a HiTrap Desalting column (GE Healthcare). The purity of the protein was analyzed through SDS-PAGE. Site-directed mutagenesis was performed using the QuikChange mutagenesis kit (Stratagene) according to the manufacturer’s instructions.

Sedimentation velocity experiments were performed on a Beckman Coulter Proteomelab XLI analytical ultracentrifuge equipped with absorbance optics. The experiments were conducted at 35000 rpm at 25 °C. Each sample (400 µl) was diluted to 0.8 AU at 280 nm in 1.2 cm optical path cell in 10 mM Hepes, 150 mM NaCl, pH 7.4. Radial absorbance scans were collected at 280 nm at a spacing of 0.003 cm in a continuous scan mode (50 scans/cell) and result was obtained over three technical replicates. Sedimentation coefficients were determined using the program Sedfit [22] (provided by Dr P. Schuck, National Institute of Health) and were reduced to water (S20,W) using standard procedures.

Equilibrium unfolding experiment of wild-type (WT) ZO1-PDZ2 was performed on a JASCO circular dichroism (CD) spectropolarimeter (JASCO, Inc., Easton, MD), in a 1 mm quartz cuvette at 230 nm. Protein concentration was typically 10 µM. The buffer used was 10 mM Hepes, 150 mM NaCl, pH 7.4 and the temperature was 25 °C. The denaturing agent was guanidine-HCl. Equilibrium unfolding experiment of PDZ2-ZO1 F207W was performed on a Fluoromax single photon counting spectrofluorometer (Jobin-Yvon, Edison, NJ, USA). The protein was excited at 280 nm, and emission spectra were recorded between 300 and 400 nm at increasing guanidine-HCl concentrations. Experiments were performed with the protein at a constant concentration of 2 μM, in buffer 10 mM Hepes, 150 mM NaCl, pH 7.4 at 25 °C, using a quartz cuvette with a path length of 1 cm. Data were fitted by a sigmoidal equation: Kinetic (un)folding experiments were performed on an Applied Photophysics SX-18 stopped-flow apparatus, monitoring the change of fluorescence emission, exciting the sample at 280 nm, and recording the fluorescence emission using a 320 nm cutoff glass filter. The experiments were performed at 25 °C using guanidine-HCl as denaturant agent.

The interaction between the C-terminal portion of the E protein (ligand) and the purified ZO1-PDZ2 domain (analyte) was measured by SPR using a Biacore X100 instrument (Biacore, Uppsala, Sweden). The N-terminal biotinylated tetradecapeptide of SARS-CoV-2 (VKNLNSSRVPDLLV) was obtained from GenScript (Piscataway, NJ, USA) and was immobilized on a Sensor Chip SA, precoated with streptavidin from Biacore AB (Uppsala, Sweden). The capturing procedure on the biosensor surface was performed according to manufacturer’s instructions setting the aim for ligand immobilization to 1000 response unit. The running buffer was Hepes-buffered saline-EP 1X, which contains 10 mM Hepes, 150 mM NaCl, 3 mM EDTA and 0.05% v/v Surfactant P20 (Biacore AB), pH 7.4. Protein domains were dissolved in the running buffer and binding experiments were performed at 25 °C with a flow rate of 30 µl/min. The association phase between ligand and analyte was followed for 180 s, while the dissociation phase was followed for 300 s. The highest concentration obtained for ZO1-PDZ2 domain was 50 μM. Concentrations measured in the SPR assay were obtained by successive dilutions: 50 μM, 25 μM, 10 μM, 5 μM, 2.5 μM, 1.25 μM, 0.5 μM and 0.25 μM. The 0.5 μM concentration was repeated twice to test the reproducibility of data. To regenerate the surface of the chip, complete dissociation of the active complex was achieved by the addition of 2 M NaCl for 30 s before the start of each new cycle. When experimental data met quality criteria, data were analyzed using Biacore X100 Evaluation Software. An affinity steady-state model was applied to fit the data, as kinetic parameters were out of the range measured by the instrument but an equilibrium signal of interaction was clearly detected. Therefore, the specific KD was determined with a confidence interval associated with a standard error value to avoid any bias.

All the structures used in this study, namely ZO1-PDZ2 monomer folded state (from AlphaFold [23]), ZO1-PDZ2 dimer (PDB code: 2RCZ [20]), and ZO1-PDZ2 dimer in complex with connexin 43 (Cx43, PDB code: 3CYY [24]), were prepared using Protein Preparation Wizard [25] from Schrödinger software, hence residues sidechain and loops have been filled in the structures when missing, then, hydrogen bond orientations have been optimized and whole structure minimized. To build the complex between SARS-CoV-2 E-SLiM and monomeric ZO1-PDZ2, the receptor grid was created by selecting βB residues with a size of 20 Å, in order to make the 14mer peptide to fit in, then the viral C-terminal tetradecapeptide, built using Maestro tools, was docked using Glide [26]. For SARS-CoV-2 E-SLiM – ZO1-PDZ2 in its dimeric state, the receptor grid was created by selecting the endogenous ligand, with the same size of the previous one, present in the structure, thus the viral peptide was docked. For both systems, the binding pose with the lowest docking score were selected and are shown in Fig. S2.

The coordinates of ZO1-PDZ2 and ZO1-PDZ2:SARS-CoV-2 E-SLiM peptide complexes were submitted to the CHARMM-GUI server [27], using the Solution Builder tool. The protein structures were centered in a cubic simulation box with a minimum distance of 10 Å between the protein and the box edge. The protein topology was generated using CHARMM36m [28] force field and the box was filled with a CHARMM-modified TIP3P water molecules [29], [30], [31]. The protein termini were treated as ionized using standard N- and C- terminal patches. All crystallographic waters were kept in the system. Next, the system net charge was neutralized by adding necessary counterions and a total ionic strength of 150 mM of KCl to mimic cellular conditions. For polarizable simulations, each system was equilibrated using the CHARMM36m force field for 1 ns using the OpenMM [32] simulation engine (version 7.6) while the coordinates non-hydrogen atoms of the protein-ligand complexes were restrained with a force constant of 500 kJ/(mol nm2). The Langevin integrator with a friction coefficient of 4 ps−1 was used with a 2-fs integration step and a temperature of 298 K, while a Monte Carlo barostat algorithm was applied to maintain the pressure at 1 bar. Short-range van der Waals forces were switched smoothly to zero from 10 to 12 Å and electrostatic interactions were calculated via the particle mesh Ewald method [33], [34] with a real-space cutoff of 12 Å. The equilibrated coordinates were converted to the Drude-2019 polarizable force field [35], [36] using the CHARMM software [37] by adding Drude oscillators and lone pairs while also converting the TIP3P water model to the polarizable SWM4-NDP model [38]. The positions of Drude oscillators were relaxed via 1000 steps of steepest descent algorithm, followed by an additional 500 steps of adopted-basis Newton-Raphson minimization. The system was then equilibrated in OpenMM for 1 ns using a Langevin integration algorithm. Temperature was regulated using a dual Langevin thermostat, friction coefficients of 20 ps−1, coupling real atoms to a thermostat at 298 K and the Drude oscillators to a low-temperature relative thermostat at 1 K. Pressure was maintained at 1 bar via the Monte Carlo rescaling algorithm. Nonbonded interaction settings were the same as described above, except that the van der Waals potential was switched to zero from 10 to 12 Å, per the convention used in the force field. The integration step was set to 1 fs and polarization catastrophe was avoided by applying the Drude “hard wall” constraint [39], allowing a maximum Drude-atom bond length of 0.2 Å. After equilibration, restraints were removed and three 500-ns production simulations were performed for further analysis. Coordinates were saved every 10 ps. Table 1. Following the number of particles present in each system: Percentages of formation of hydrogen bonds interactions along MD trajectories were calculated by custom procedures using the MDAnalysis Python library [40]. Polar and hydrophobic interactions between peptides and receptor residues were modelled as average coordination numbers via a continuous, differentiable switching function: With the I and j indexes running over the interacting atoms within the peptide fragment of interest and those within a chosen receptor residue, with 0 ≤ strength ≤ ni nj (where ni and nj are the total number of atoms of the peptide fragment and the receptor residue, respectively, able to make the chosen type of interaction). For hydrophobic interactions, only carbon atoms are considered, for polar interactions only oxygen and nitrogen atoms. The same “strength” function was calculated for salt bridge, cation-π and π-π interactions. In this case ni = nj = 1 because virtual atoms are defined at the center of mass of aromatic rings or charged groups, hence 0 ≤ strength ≤ 1. For hydrophobic, salt bridge, cation-π, and π-π interactions a = 6, b = 12 and r0 = 6 Å, 5 Å, 4 Å, and 5.5 Å respectively, while for polar interactions a = 8, b = 12, r0 = 2.5. The chosen values of r0 account for the typical interaction distance plus an amplitude associated with thermal motion (e.g. ∼4.5 Å + ∼1.5 Å in the case of hydrophobic interactions) [41]. These coordination numbers (one per residue of the binding site) can be calculated on single structures as well as averaged along trajectories. Dipole moments of the various residues and functional groups in the protein structures were calculated using built-in CHARMM functions [37]. RMSD-based clustering was performed using CHARMM software setting the radius parameter to 1.2 nm and analyzing the last 250 ns of each simulations. As regard CHARMM36m E-SLiM in complex with ZO1-PDZ2 in its monomeric state, just the 2 simulations were the detachment doesn’t occur have been analyzed.

To determine whether ZO1-PDZ2 was monomeric or dimeric under our experimental conditions, we performed analytical ultracentrifugation analysis. Results are reported in Fig. 1D. As described in the Introduction, a complex behavior regulating ZO1 binding has been described [42], in which ZO1 undergoes dimerization through domain swapping occurring at the level of the PDZ2 domain, with the consequent formation of two distinct, PDZ-shaped, fully functional binding modules. The experiment was performed on an engineered pseudo wild-type ZO1-PDZ2, namely F207W, which was then employed in folding experiments (see below in the text). Our data clearly show only one symmetric peak (97%) with a sedimentation coefficient of S20,w of 1.9 S at 25 °C, in accordance with the molecular weight of the monomeric form of ZO1-PDZ2 (∼10 kDa), unequivocally showing that it is a monomer in the experimental conditions used. We further analyzed the folding properties of ZO1-PDZ2 via equilibrium and kinetic (un)folding experiments. Equilibrium (un)folding experiments were conducted by challenging a constant concentration of the ZO1-PDZ2 F207W with increasing concentration of denaturant (guanidine-HCl) and following the change of the intrinsic fluorescence as function of denaturant concentration (Fig. 1 A). The buffer used was 10 mM Hepes, 150 mM NaCl, pH 7.4 and the temperature was set to 25 °C. A single, sigmoidal transition between native and denatured states was observed. Notably, the mD-N value (an index of the cooperativity of the reaction, defined as δΔG/δ[denaturant] and correlated with the change of the accessible surface area upon unfolding) was 2.2 ± 0.2 kcal mol−1 M−1, as expected for a protein of 92 residues [43] with a ∆GD-N = 4.4 ± 0.3 kcal mol−1. To check the effect of the F207W mutation on the thermodynamic stability of the protein, we conducted equilibrium unfolding on the wild-type variant by following CD signal at 230 nm at different concentrations of denaturant. The denaturation curve is reported in Fig. 1 B. The result of fitting was mD-N = 2.0 ± 0.2 kcal mol−1 M−1, and ∆GD-N = 3.6 ± 0.4 kcal mol−1. These thermodynamic parameters appear comparable to the ones obtained by fluorescence, concluding that the F207W mutation does not influence the stability of the PDZ2 domain. To increase our understanding of the folding properties of ZO1-PDZ2 domain, we performed folding and unfolding kinetic experiments through stopped-flow methodology (SX-18 Applied Photophysics) under the same experimental conditions used for equilibrium experiments. Importantly, at all denaturant concentrations, folding and unfolding time courses were satisfactorily fitted with a single exponential equation; typical unfolding and refolding traces are reported in Fig. 1E,F. The logarithm of observed rate constants (k obs) were plotted versus denaturant concentrations (chevron plots). Fig. 1 C highlights the clear presence of a “roll-over effect” in the unfolding arm of the chevron plot (i.e., a deviation from linearity in the dependence of the logarithm of the k obs). This observation is classically interpreted as a signature of complex folding, reflecting the presence of a folding intermediate along the reaction pathway, either a high-energy intermediate, implying a change in the rate-limiting step at low and high denaturant concentrations, or a low-energy accumulating intermediate [44], [45]. Together, these results indicate that the PDZ2 domain is capable of folding in its monomeric state.

Numerous previous studies have concluded that coronavirus E proteins interact with different PDZ domains [7], [8], [46], [47]. The interaction between ZO1-PDZ2 domain and WT SARS-CoV-2 E-protein was previously described by Shepley-McTaggart et al. [21]. These authors reported a KD value of 29 nM by performing Homogenous Time Resolved Fluorescence (HTRF). This reported affinity value is quite high compared with the affinity existing between antigen and antibody [48] and also higher than the affinity measured for SARS-CoV-2 E-protein and another member of the PDZ family, PALS1 PDZ8). Based on the findings presented, we decided to measure the interaction between ZO1-PDZ2 domain and WT SARS-CoV-2 E protein using a different experimental technique, namely SPR. To achieve this goal, N-terminal biotinylated tetradecapeptides of the C-terminal E protein were immobilized on a streptavidin chip (SA) and used as ligands, while the ZO1-PDZ2 domain was injected as the analyte in our SPR assay in a wide range of final concentration [0.5 – 50 µM] using a multi-kinetic mode. We found that ZO1-PDZ2 domain binds the C-terminal sequence of SARS-CoV-2 E protein with a KD = 2.1 ± 1.1 μM, a much lower affinity than that reported by Shepley-McTaggart and coworkers [21]. The SPR results ( Fig. 2) confirm that PDZ domains are suitable partners for coronavirus E-PBM, even if our result reports an equilibrium dissociation constant in the order of micromolar range. Despite conflicting results, the KD value obtained in this study seems more plausible. In fact, binding affinity for ZO1-PDZ2 domain with the SARS-CoV-2 E protein is in good agreement with other measurements involving SLiMs and PDZ domains, as previously mentioned [8]. Moreover, a similar KD value for this analyte with its endogenous ligand, Cx43, was also determined [24].

Since our recombinant protein was characterized and shown to be a monomer, a computational model for this conformation was built, assuming that ZO1-PDZ2 may exist in two different conformational states: a canonical folded state and an extended, partly unfolded state, that will lead to the formation of the dimer as reported in the PDB (entries 2RCZ and 3 CYY) [20], [24]. The predicted structure of the human ZO1 protein was retrieved from the AlphaFold Database [23], [49] and residues 184–262, comprising the PDZ2 domain, were extracted to be compared with the ZO1-PDZ2 domain present in the crystal structure of the dimer (PDB 2RCZ20). The ZO1-PDZ2 monomer structure predicted by AlphaFold scores well in terms of both lDDT (local Distance Difference Test) analysis and Predicted Alignment Error (Fig. S1 A, B). Nearly all of the residues have an lDDT score above 80 and the relative distances between residues are all well resolved. The structure was further analyzed with ProSA-web [50] tools and its results affirm the strength of the model predicted by Alphafold, obtaining a z-score of − 4.9 (Fig. S1 C) and a satisfying Local Model Quality plot (Fig. S1 D) where the energies remain negative throughout the sequence in the 40-residue average. In the predicted model of the monomer, the sequence involved in the domain swapping (residues 186–206) is folded back onto the PDZ core ( Fig. 3C). Specifically, residues 200–203 form the βB strand (cyan), antiparallel to βC, while residues 191–199 form the connecting loop (green) to the residues 186–198, which are now arranged in the βA strand (blue) antiparallel to βF. To assess the stability of the ZO1-PDZ2 in the two conformations, monomer (Fig. 3C) and dimer (Fig. 3A), three 500-ns MD simulations were performed with the Drude-2019 force field. The structural parameters derived by the analysis, i.e. RMSD, α-helix and β-sheet contents were obtained considering the monomeric state in comparison with a single polypeptide chain of the dimer of ZO1-PDZ2. The Cα RMSD of the ZO1-PDZ2 monomer plateaued in the second half of each of the 500-ns simulations in all three replicas after an initial increase. In contrast, the dimer remained stable in the entire time window maintaining its initial structure as derived from X-ray crystallography (PDB code: 2RCZ20) ( Fig. 4 A,B). α-helix and β-sheet content time series, obtained via DSSP analysis [51], are shown in panels C-F of Fig. 4. The ZO1-PDZ2, both in its monomeric and dimeric forms, had the constant α-helical content over the course of each simulation. The β-sheet content of dimeric ZO1-PDZ2 was constant over time, whereas in the monomeric form, there was a small reconfiguration and organization of novel intrachain β-strand units, leading to a slight increase in β-sheet content (Fig. 4E). Indeed, the predicted structure of the ZO1-PDZ2 monomer contains novel folded β-strands (βA and βB) that would have to disorder and elongate to engage in the domain swapping observed in the experimentally resolved dimer structure. These β-strands interact with the PDZ2 core domain via an intrachain hydrogen bonding network that is shown in Fig. 5A. Interestingly, this hydrogen bonding network involves the same pattern of interactions as in the dimer, since the interacting residues are from the βA and βB strands of the other monomer (Fig. 5B). Thus, the two end states may have very similar energies, but additional studies would be necessary to test this hypothesis and identify any free energy barriers between the monomeric and domain-swapped states. To compare the flexibility of monomeric and dimeric ZO1-PDZ2, we computed the root-mean-square fluctuation (RMSF) for Cα atoms ( Fig. 6). The flexibility of the ZO1-PDZ2 core (residues 210–262) was not significantly affected by the three-dimensional arrangement, with similar RMSF values in both states, whereas residues 184–210 behaved differently. Novel βA strand (residues 185–190) in the monomeric folded conformation had lower RMSF values compared to the same sequence in the dimer, despite being stabilized by a very similar set of hydrogen bonds in both systems (Fig. 5). Both in the monomeric and dimeric forms, this initial sequence is followed by a more flexible region (residues 191–196) that corresponds to the loop connecting βA to βB in the monomer, or to the loop participating to the domain swapping domain in the dimer (see green regions in Fig. 3B,C). In the sequence spanning βB to βC in the PDZ core domain (residues 197–209), the monomer was always more rigid than the dimer. Even the short loop connecting βB to βC (residues 202–206) was not as flexibile as the corresponding residues in the dimer. As the process of domain swapping involves changing the microenvironment around βA, βB, and βC structural motifs, we also investigated the peptide-bond and sidechain dipole moments in these systems. The explicitly polarizable nature of the Drude force field allows for the electronic plasticity of the structures to be modeled, and we sought to determine if the electronic environment in the monomer and dimer of ZO1-PDZ2 differed. Differences in peptide-bond and sidechain dipole moments are shown in Fig. 7. Few differences between monomer and dimer conformations were observed, though they generally manifested in the region encompassing βB to βC. As expected, in the monomer, the groove left by the domain swapping is perfectly replaced by its own refolded strand, managing to combine perfectly with the PDZ core. Interestingly, we observed a pronounced polarization effect in the short loop corresponding to the flexible region of residues 202–206. Peptide bonds in this loop were more polarized in the monomer (Fig. 7A), due to the backbone being completely exposed to the aqueous solvent, whereas in the dimer conformation this same region is buried due to the domain swapping. Correspondingly, the dipole moments of the constituent sidechains were strongly depolarized in the monomeric state (Fig. 7B), suggesting a positioning of these sidechains toward the inner core of the protein domain. On the contrary, the sidechains of the dimer for the same region result more polarized as they were completely exposed to the aqueous solvent.

Our experimental and computational data clearly show that the ZO1-PDZ2 recombinant protein is also stable in a monomeric state. Moreover, the SPR assay also demonstrated that the monomer is still able to interact with SARS-CoV-2 E-SLiM. Accordingly, we focused our attention on identifying structural determinants that characterize this interaction and we subsequently simulated this heterodimeric structure. The E-SLiM tetradecapeptide was docked onto the monomeric ZO1-PDZ2 domain in its putative binding pocket between βB and αB (see Fig. S2 A). Successively, the resulting complex was challenged via three 500-ns MD simulations using the Drude-2019 force field. The frequencies of different interaction types among key residues are shown in Fig. 8 A, and the most representative snapshot is shown in Fig. 8 panels B and C. The SARS-CoV-2 E-SLiM peptide largely interacted with the ZO1-PDZ2 domain through a combination of favorable interactions. The negatively charged carboxylate group of the C-terminal V75 residue is involved in a complex network of hydrogen bonds with the protein backbone moieties of Y198, G199, and L200. This network is necessary to orient in a β-sheet-like arrangement adopted by the D72, L73, and L74 of the E-SLiM, stabilizing this antiparallel-like conformation through a hydrogen bonding network with the protein partners R201 and L202 (Fig. 8 A, C). This backward positioning of the SLiM peptide, mainly caused by the hindrance of the Y198 ZO1-PDZ2 sidechain, allow the V75 sidechain to reside deeply in the ZO1-PDZ2 hydrophobic pocket to interact with A245 and I249. Interestingly, viral D72 binds via salt bridge interaction both human protein residues R201 and K209, meanwhile K246 remains close to P71 carbonyl oxygen due to polar interaction. We compared these results with a second set of simulations performed with the same engine (OpenMM32) but different force field (CHARMM36m28) to determine if there was any difference between polarizable (Drude) and additive (CHARMM36m) force fields. We again performed three 500-ns simulations. In one of them, the E-SLiM peptide partially dissociated from the ZO1-PDZ2 binding partner ( Fig. 9B) and in the other two simulations the complex between viral peptide and PDZ domain was present during the entire simulation time even if it was less stable, as quantified in terms of Cα-RMSD (Fig. 9 A, B). In fact, the E-SLiM peptide moved upward and outward during the entire simulation time, with C-terminus of V75 trying to interact with ZO1-PDZ2 residue K191 to form a salt bridge (Fig. 9C, green structure). Thus, in these simulations, V75 preferred this single ionic interaction over a more stable hydrogen bonding network with the triad of PDZ residues Y198, G199 and L200 of ZO1-PDZ2. Indeed, in the Drude simulations, the E-SLiM remained very close to its starting configuration (Fig. 9A), thanks to this hydrogen bond network with the backbone of ZO1-PDZ2 βB. It has previously been observed that many additive force fields overstate the strength of salt bridge interactions [52], and in a complex interface such as the one that arises between ZO1-PDZ2 and E-SLiM, this property may be more accurately modeled with a polarizable force field like Drude-2019. Moreover, the last 4 residues form a structured β-strand-like conformation for the 71% of the simulation time. This state is never sampled in CHARMM36m simulations, where the peptide remains fully unstructured along the entire simulation. In Fig. 9C,D, the two most representative conformations for both set of simulations are shown for comparison. To further investigate the role of the force field in the formation of a protein peptide complex, we also prepared another set of simulations to investigate the binding mode of the E-SLiM peptide with dimeric ZO1-PDZ2. To build this complex, the viral peptide was docked into both Cx43 (the endogenous peptide partner) binding sites in ZO1-PDZ2, starting from the three-dimensional structure reported in PDB 3CYY [24]. Three 500-ns MD simulations of these complexes were performed using both the Drude-2019 and CHARMM36m force fields (starting conformation in Fig. S2 B). In the Drude simulations, the last 4 C-terminal residues of the E-SLiM peptide interacted with the ZO1-PDZ2 dimer in a similar fashion as in the monomer, where the βB strand present in the monomeric conformation is replaced by the strand of the other monomer performing domain swapping. In fact, the C-terminal carboxylate group of V75 was always involved in the hydrogen bonding network with the same residues Y198, G199 and L200 of the PDZ domain. Further, the following 3 residues of the viral peptide are organized in β-sheet-like conformation, which allows the E-SLiM R69 to engage a complex salt-bridge network (see Fig. 10). Moreover, the proline residue at position 71 forces the peptide to turn, facilitating the formation of an intramolecular salt bridge between D72 and R69 and an intermolecular salt bridge between R69 and the ZO1-PDZ domain residue E238. Finally, D72 interacted with K209 and R201, which in turn was involved in a salt bridge with E210 and S68 also stabilized the binding pose via polar contact with R201. This binding pose of the SARS-CoV-2 E-SLiM is very similar to that of Cx43. In fact, Cx43 and E-SLiM have similar sequences, as shown in Fig. 11 where the endogenous peptide shows a proline residue after the PDZ binding motif, which is followed by an arginine and preceded by an aspartate. This motif can be also found in SARS-CoV-2 E-SLiM peptide with the only difference being a valine residue between the P71 and R69. In the CHARMM36m simulations, the viral peptide C-terminal portion, even in the binding site of the dimeric conformation of the ZO1-PDZ2 domain, never adopts a β-strand-like conformation. As a consequence, a main difference between Drude and CHARMM36m simulations is the number of polar contacts observed, which were fewer using the CHARMM36m force field as demonstrated by the comparison of the interaction histograms ( Fig. 12A,B). Moreover, the E-SLiM C-terminal V75 carboxylate group did not interact with Y198, G199, and L200 via hydrogen bonds as in Drude simulations, seeking to form salt-bridge interaction with available partners such as K191 and K246 (Fig. 12D). However, the E-SLiM peptide never shifts further reaching K191, as previously seen in the interaction with the monomeric conformation of the human PDZ domain, due to the salt bridge between E238 and R69 that holds back the viral peptide.

The reason of the success of the virus causing Covid-19 pandemic to spread worldwide has been largely debated [53], [54], [55]. Among possible causes of its ability to easily diffuse in the human organism, a broad consensus exists for the presence of a highly conserved C-terminal portion of the E protein containing a PDZ-binding motif, which allows it to interact with host cell proteins with one or more PDZ domains [7], [8], [56]. In this study, we focused our attention on one of the main targets among these PDZ domains, namely ZO1-PDZ2, as it has been previously reported a very high affinity (in the nanomolar range21) between the human protein and the viral E-SLiM. A first clear result of our study is that this ZO1-PDZ2 domain may exist as a stable monomer, that is alternative to the dimeric form, largely present in the human cells, where it is necessary to drive its function. Moreover, we clearly established that even this monomeric species can bind the E-SLiM, thus suggesting that it may also be targeted during viral invasion of the host cells, where the monomer may be transiently present during the cell cycle being even more easily targeted by the viral peptide. Specifically, the monomeric species that is formed in the endoplasmic reticulum could be hijacked by interacting with the newly synthesized E-protein after viral infection, thus preventing the correct trafficking of ZO1, reducing the maturation and tight-junction formation [57]. Moreover, also a possible direct mechanistic cause for the damaging of the airway epithelia, through disruption of tight junction function, may be the consequence of a competition between a functional dimeric form of ZO1-PDZ2 domain and a mis-functional monomeric species that is stabilized through the binding with the E-SLiM peptide [58]. Moreover, the computational analysis of the complex between the viral E-SLiM with ZO1-PDZ2, both in its monomeric or dimeric states, demonstrated that the binding interface between the two proteins is sensitive to electronic polarization as it involves a network of potentially competing hydrogen bonds and salt bridges. The Drude-2019 polarizable force field modeled these interactions preserving the predicted binding pose, while the CHARMM36m force field did not model these interactions as well. This effect is probably caused by a bias toward salt-bridge interaction over the oriented polar network, that is formed during the binding process. Thus, the Drude-2019 polarizable force field may be more suitable for studying non-structured peptide binding behavior, where the gain of structure upon binding is a key structural element.

Noah Giacon: Investigation, Formal analysis, Writing – original draft. Ettore Lo Cascio: Investigation, Formal analysis, Writing – original draft. Darcy S. Davidson: Formal analysis, Writing – original draft. Marcelo D. Polêto: Formal analysis, Writing – original draft. Justin A. Lemkul: Writing – review & editing. Valeria Pennacchietti: Investigation. Livia Pagano: Investigation. Carlotta Zamparelli: Investigation. Angelo Toto: Conceptualization, Writing – review & editing. Alessandro Arcovito: Conceptualization, Writing – review & editing.

On behalf of all the authors whose names are listed immediately below, I certify that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers’ bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent-licensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript.