The starting structures for MD simulations were either specific frames of deposited PDB files or partially modeled to avoid structural gaps. The structures with the highest resolution and sequence coverage available were chosen for both Rec and NCS1, resulting in the following models: Ca²⁺-free Rec (“tense Recoverin” or Rec-T, Tanaka et al., 1995), Rec with one Ca²⁺ bound to EF3 (“intermediate state Recoverin” or Rec-I, Ames et al., 2002), Ca²⁺-loaded Rec (“relaxed Recoverin” or Rec-R, Ames et al., 1997), Ca²⁺-loaded Rec bound to GRK1 peptide (Rec-GRK1, Ames et al., 2006; Zernii et al., 2011; Ames and Lim, 2012), Ca²⁺-loaded uncomplexed NCS1 (“isolated NCS1” or NCS1-iso), Ca²⁺-loaded NCS1 bound to D₂ Dopamine receptor peptides (NCS1-D₂R) and Ca²⁺-loaded NCS1 bound to Rhodopsin Kinase peptide [NCS1-GRK1, (Bourne et al., 2001; Pandalaneni et al., 2015)]. Details for molecular modeling are provided in Supplementary Methods.

MD simulations of Rec and NCS1 states were performed using GROMACS 2016.1 simulation package (Abraham et al., 2015) and CHARMM36m (Huang et al., 2017) all-atom force field, where parameters for N-myristoylated Gly were generated manually (available upon request). Simulations were performed in a dodecahedral box with Periodic Boundary Conditions applied, where proteins were located at 12 Å distance from box boundaries. Solvent was modeled as TIP3P water, system was neutralized first with 1 mM MgCl₂, then after addition of 150 mM KCl. The size for each simulated system was the following: 33,837 atoms for Rec-T, 38,878 atoms for Rec-I, 56,010 atoms for Rec-R, 86,456 atoms for Rec-GRK1, 46,809 atoms for NCS1-iso, 42,882 atoms for NCS1- D₂R and 50,584 atoms for NCS1-GRK1.

All structures underwent substantially the same pre-production steps as previous studies concerning GCAP1 (Marino et al., 2015b), briefly consisting of steepest descent (F = 1,000 kJ/mol^*nm) and conjugate gradients (F = 500 kJ/mol^*nm) energy minimization, then equilibration at 310 K for 2 ns backbone-constrained and 2 ns unrestricted MD simulations in NVT ensemble. After equilibration, 200 ns unrestrained MD simulations were performed in NPT ensemble at 310 K and 1 atm for each system. Equilibration and production MD simulations were independently replicated 5 times by changing the random seed for initial velocity generation as previously described (Marino and Dell'orco, 2016), to achieve exhaustive sampling of the conformational space.

Collective protein motions were identified using Principal Component Analysis (PCA) after diagonalization of the Cα covariance matrix (Amadei et al., 1999) calculated on single 200 ns replicas and concatenated 1 μs trajectories for each protein form, after superimposition to the final structure of the equilibration phase of replica one. Eigenvalues and eigenvectors extracted from the diagonalized matrix represent amplitude and direction of collective motions (Principal Components, or PC) and are ranked in decreasing order of eigenvalues, meaning that the first PC (PC1) accounts for the largest collective motion, PC2 for the second largest and so on.

The Essential Subspace (ES) of the first 20 PC, accounting for 78–93% of total motion of Rec and NCS1 states, was identified as in Marino and Dell'orco (2016) to compare the conformational space sampled by each replica and 1 μs trajectories using the Root Mean-Square Inner Product (RMSIP) index (Amadei et al., 1999). Details for cosine content (c₁) (Hess, 2002) and RMSIP index calculations are provided in Supplementary Methods.

The sampling convergence of MD simulations was also assessed by performing Linear Discriminant Analysis (LDA, Martinez and Kak, 2001) on the projection of the frames of the 5 replicas on PC1 and PC2 calculated on the concatenated 1 μs trajectories. This supervised classification method decreases the number of features describing data and combines them to find a hyperplane that maximizes separation between means of projected classes and minimizes the variance within each projected class. In our case a two-dimensional feature space (PC1 and PC2) was classified using a single Linear Discriminant (LD1) and if a single conformation from a replica was classified as possibly belonging to a different replica, then the two replicas could be considered as consistent and therefore concatenable. Data presented in Figures S3, S4 was smoothed using Kernel Density Estimation smoothing (Rosenblatt, 1956; Parzen, 1962).

Dynamic information from concatenated 1 μs trajectories was converted into a PSN using PyInteraph software (Tiberti et al., 2014), with the same parameters as in Marino and Dell'orco, 2016, the only exception being the mass of atoms defined according to CHARMM36m force field (Huang et al., 2017). Briefly, non-bonded interaction between side chains (H-bonds, electrostatic and hydrophobic interactions) was computed as a percentage of frames where distance and angle constraints were satisfied. Interactions were filtered using the hydrophobic cluster size criterion (Vishveshwara et al., 2009), therefore the size of the biggest hydrophobic cluster was calculated at 0.1% persistence intervals, and the persistence threshold (p_T) was computed as in Marino and Dell'orco (2016) and rounded at the lowest decimal value. For each case p_T was calculated and all interactions above the threshold were joined in the unweighted PSN graph representing the specific protein state.

The degree of connectivity was calculated for all residues in each PSN, amino acids with ≥ 4 connections were considered hubs and reported in Tables S1–S3. The Communication Robustness index (CR, Marino and Dell'orco, 2016) of two residues permits the evaluation of the structural information transfer between amino acids within a protein three-dimensional structure, and to highlight potential allosteric mechanisms. Details for CR calculation are provided in Supplementary Methods. CR index was therefore computed to measure the communication between Rec bidentate Ca²⁺-coordinating residues E85 and E121 and GRK1 interface residues (Table S4, Ames et al., 2006; Zernii et al., 2011). Analogously for NCS1, CR index was calculated between bidentate Ca²⁺-coordinating residues E84, E120, and E168 and with GRK1 or D₂R interface residues (Table S4, Pandalaneni et al., 2015). For the sake of clarity, only interface residues with the 5 highest CR values for each presented state were reported in Figures S5, S6. Moreover, intramolecular communication was evaluated by analyzing CR calculated between the myristoyl group of Rec and GRK1-interface residues (Figure S7A) and among the previously mentioned Ca²⁺-coordinating residues of Rec and NCS1 (Figure S7B), to investigate intramolecular communications between EF-hand motifs.

Finally, also intermolecular communication was investigated by computing CR index between GRK1 peptide residues that were structurally solved (L6-I16) complexed with both Rec and NCS1 and their respective EF-hand representatives (Figure S8).

For visualization purposes, pathways connecting EF-hand representatives and GRK1 residues with the highest CR were chosen according to the highest cumulative Selective Betweenness (SB, Marino and Dell'orco, 2016). Details for SB calculation are provided in Supplementary Methods.

Protein sequences for bovine Rec (Uniprot: P21457), human NCS1 (Uniprot: P62166), human GCAP1 (Uniprot: P43080) were subjected to Multiple Sequence Alignment (MSA) with Clustal Omega (Sievers et al., 2011) to evaluate residue conservation among homologous NCS. Moreover, each protein sequence was also subjected to MSA with Clustal Omega with up to 250 Uniprot sequences having at least 50% sequence identity (s.i.) and 80% sequence coverage with respect to the longest sequence in the cluster. In detail, Rec was aligned with 122 sequences from Uniref50_P21457 and NCS1 was aligned with 250 sequences from Uniref50_P62166. Human GCAP1 was instead aligned with 101 sequences from Uniref50_P46065, referring to the bovine sequence, constituting the seed for the linked UniRef50 cluster (s.i. 94%).

Residue conservation in UniRef50 clusters shown in Figure 6 and in Tables S1–S3, S5 was calculated as the ratio between the number of sequences where the residue of the seed sequence was conserved and the total number of sequences in the cluster. When a residue was a hub in different signaling states, even with a different degree, only one occurrence of the residue was considered for the calculation of the average of a given hub degree. Average values for each protein were mediated and subjected to linear regression after passing Shapiro-Wilk (Shapiro and Wilk, 1965) (p = 0.05) and constant variance (p = 0.05) tests, with R² = 0.62.

Hub degree shown in Figure 5 represents the highest hub degree exhibited by a residue in one of the simulated states for each protein, data about GCAP1 are taken from Marino and Dell'orco (2016).

Pairwise sequence alignments were performed with Needleman-Wunsch (Needleman and Wunsch, 1970) algorithm.

Rec-T, Rec-I, Rec-R, and Rec-GRK1 were simulated to evaluate the network of persistent interactions and any potential difference in the inter- and intra-molecular allosteric properties. Such differences are ascribable to the physiological conformational transition of the sensor as a consequence of the decrease of intracellular [Ca²⁺] upon activation of the phototransduction cascade (Ames et al., 1996). On the other hand, NCS1-iso, NCS1-GRK1, and NCS1-D₂R were simulated to evaluate differential communication of the same protein in the absence and in the presence of different targets. Furthermore, the presence of the same target, namely GRK1, interacting with both Rec and NCS1 allowed for a comparative analysis of the communication routes between Ca²⁺-binding sites and a common target.

For each protein state 5 × 200 ns replicas of MD simulations were run; individual replicas and 1 μs concatenated trajectories were subjected to PCA of Cα motions to assess consistency and exhaustiveness, as previously described (Marino and Dell'orco, 2016). Cosine content c₁ of PC1 was the first necessary condition to achieve sampling convergence, as it was demonstrated that values close to 1 are indicative of insufficient sampling (Hess, 2000; Papaleo et al., 2009). The c₁ for the concatenated trajectories of Rec ranged between 0.0006 and 0.1244, those for NCS1 ranged between 0.0150 and 0.1597, while the average c₁ of the five replicas ranged between 0.4670 and 0.5792 for Rec and between 0.1421 and 0.3302 for NCS1. The decrease of c₁ values from the average of the single replicas to that of the concatenated trajectories suggests that the longer timescale allows for a more significant and exhaustive conformational sampling.

The RMSIP of the first 20 PC of the ES highlighted a substantial overlapping, and therefore reproducibility, of the conformational subspaces sampled by each individual replica and by the concatenated 1 μs trajectories, as shown by RMSIP values for NCS1 and Rec, ranging from 0.763 to 0.895 (Figure S1), and from 0.665 and 0.828, respectively (Figure S2).

LDA was performed on the projection of MD trajectories along PC1 and PC2 for each protein state. Results showed that for each couple of replicas in all NCS1 (Figure S3) and Rec (Figure S4) states there was at least one conformation that could belong to two different replicas, thus allowing these replicas to be treated as one concatenated trajectory.

Overall, all results confirm that individual 200 ns MD replicas are independent and consistent representations of the dynamic structural behavior of both Rec and NCS1 in their different signaling states. Moreover, concatenated 1 μs trajectories exhaustively sample each protein conformational space and can be subjected to further detailed PSN analyses.

Non-covalent interactions are crucial for defining protein three-dimensional structure, yet hydrophobic, electrostatic, and H-bond interactions can have a transient or more persistent essence, depending on whether they are involved in structural information transfer or not. Therefore, we monitored the persistence of these interactions throughout the concatenated 1 μs trajectories of Rec and NCS1 signaling states and determined their significance in shaping the specific PSN (Vendruscolo et al., 2002; Atilgan et al., 2004; Papaleo, 2015) for each protein state. In the resulting graph residues are nodes and persistent interactions represent edges.

Dynamic information from 1 μs MD simulations was encoded in a static PSN (see Methods) and the analysis of graphs identified the presence of key residues (hubs) persistently interacting with many close-contact amino acids (Tables S1–S3). To discriminate whether the identified interactions were significant or transient, p_T was calculated for each trajectory to filter out short-lived interactions (see Methods) and was found to range between 15.7 and 17.6 for Rec states and between 36.5 and 37.2 for NCS1 variants.

To assess which ones among the interface residues are responsible for transmitting structural information regarding the Ca²⁺-loading state to their target, the CR index was calculated for each residue of NCS1 (Figure S5) and Rec (Figure S6) target interfaces with each respective Ca²⁺-coordinating Glu residues chosen as representative for each EF-hand motif.

Results for NCS1, summarized in Figure 1, highlight that in the absence of target peptides (orange arrows) there is a long-range communication from EF2 (Figures 1A,D) to both EF-hands of the C-domain. In detail, L101 (Figure 1A) and Q130 (Figure 1D) belong to the entering (hE) and exiting (hF) helix of EF3, respectively, and F169 (Figures 1A,D) is located at helix hF of EF4 (hF4). On the other hand, communication from EF3 (Figures 1B,E) and EF4 (Figures 1C,F) is overall short-ranged in the absence of the target, reaching for both interfaces only residues located in adjacent EF2 or EF4. Specifically, E120 communicates with F64 and V68 (Figure 1B) or with V68 and F72 (Figure 1E), which are all located on helix hE2. The EF4 representative residue E168 exhibits robust communication with M121 and I124 of helix hF3 (Figure 1C) or W103 and Q130, located, respectively, on helix hE3 and hF3.

In the presence of GRK1 (blue arrows) all three EF-hands show allosteric intradomain properties: EF2 shows robust communication with F55 of helix hE1 and F169 of helix hF4. (Figure 1A); EF3 exhibits high CR values with M131 of helix hF3 and to a lower extent with F55 of helix hE1; finally, EF4 robustly communicates with L107, located on the adjacent helix hE3, and with Y52 of helix hF1.

The pattern of structural information transfer appears different in the presence of another biological target. Overall, communication of NCS1 EF-hands with D₂R interface travels shorter distances (Figures 1D–F, blue arrows). In detail: EF2 communicates with hE1 residue F56 and hE3 residue Q130 (Figure 1D); EF4 communicates with hE3 residue L107 and hF3 residue I128 (Figure 1F); finally, EF3 communication reaches Y129, also belonging to hF3 helix, and the distant W30 of helix hE1, the only residue exhibiting long-range communication (Figure 1E). In summary, allosteric communication between EF-hands and target interface is a peculiar feature of NCS1-GRK1 complex and of the EF2 motif of NCS1-iso, while robust communication is on average shown on a mid-to-short range by NCS1-iso EF3 and EF4 and by NCS1- D₂R complex.

Results for Rec, summarized in Figure 2, highlighted robust short-range communication from EF2 to GRK1 interface in the absence of target. In Rec-T (Figure 2A) the absence of Ca²⁺ ions bound to the EF-hands is communicated from EF2 to F23 (hE1) and Y53 (hF1). In Rec-I (Figure 2B), where only EF3 is occupied by a Ca²⁺ ion, the communication of unoccupied EF2 is robust with F57 (hF1) and L90 (hF2). Finally, in Rec-R (Figure 2C), with both EF-hands occupied by Ca²⁺, again EF2 communicates with hE1 residue F23 and hF2 residue L90.

The presence of the GRK1 target changed somewhat the scenario. Rec-GRK1 complex, indeed, revealed a strikingly high allosteric communication between EF2 (E85) and hF4 residues E189 and K192 (Figure 2D), a residue previously identified as crucial for the interaction with GRK residue F3 (Zernii et al., 2011). Long-range communication was observed also when originated in EF3, as shown by the high CR exhibited by Rec-T L90, located in hF2 (Figure 2A), the previously mentioned Y53 in Rec-I (Figure 2B), hE1 residue F35 in Rec-R (Figure 2C), and hF1 residue F49 in Rec-GRK1 complex (Figure 2D). Surprisingly, regardless of the presence of Ca²⁺ or GRK1, EF3 showed a significant robust communication with the first residue of helix hF2, namely Y86. Moreover, Y86 counterpart in NCS1 (F85) was identified as one of the most robustly communicating interface residues with EF2 and EF4, particularly in the absence of targets (Figure S5).

The N-terminal myristoyl group of Rec was also found to robustly communicate with GRK1 interface residues (Figure S7A) in the absence of Ca²⁺ (Rec-T), but this robustness decreased on average upon Ca²⁺-binding to EF3 (Rec-I). Moreover, no CR index could be calculated upon Ca²⁺-binding to EF2 (Rec-R and Rec-GRK1) due to the extrusion of the fatty acid from the hydrophobic crevice via myristoyl switch mechanism (Tanaka et al., 1995), which prevented any persistent interaction with the protein.

Finally, CR index was calculated among EF-hands of Rec and NCS1 (Figure S7B), highlighting substantial medium-range information transfer between EF2 and EF3 of Rec only in the absence of Ca²⁺. NCS1, on the other hand, showed no significant communication between adjacent EF-hands, but the NCS1-GRK1 complex exhibited a surprisingly high robustness in the allosteric communication between EF2 and EF4 (Figure S7B).

Intermolecular communication between GRK1 and EF-hands of Rec and NCS1 was assessed via CR index to compare how the same target could be regulated by two proteins with overlapping interfaces. Results (Figure S8) highlighted that GRK1 residue L6 communicated specifically with NCS1-EF4, T8 with NCS1-EF3, A11 with Rec-EF3 and S13 with NCS1-EF2. Interestingly, residue N12 of GRK1 was found to communicate with both Rec-EF2 and NCS1-EF3, while F15 with EF2 in NCS1 and EF3 in Rec. Overall, the analysis highlighted that the robustness of NCS-target interaction is highly protein-dependent in terms of specificity of the signaling state.

In addition, the previously mentioned GRK1 residues with highest CR were investigated from a topological standpoint concerning the most probable PSN pathways that would allow information transfer to/from EF-hands. The most probable pathways were chosen according to the highest collective SB, nevertheless in some cases multiple pathways with the same score were reported (Figure 3). Pathway analysis (Figure 3) indicated that the pathway from NCS1 EF2 to S13 (Figure 3A; Video V1) included Ca²⁺ ion and that the three GRK1 residues of the pathway, namely L6, V9 and S13, were also involved in the communication of the peptide with NCS1 EF3 and EF4. The pathway from NCS1 EF3 to N12 was the most branched out (Figure 3B; Video V2) with two forks and a double intermolecular communication required to correctly route Ca²⁺-binding information to V9. In addition, EF3-N12 pathway shared the GRK1 residues V9 and S13 with EF2-S13 pathway and both GRK1-I16 and NCS1-Y108 residues with I16 EF4-L6 pathway (Figure 3C; Video V3). Pathway EF4-L6 shared GRK1 residue L6 with EF2-S13 pathway.

The same pathway analysis performed on Rec EF2 with GRK1-N12 (Figure 3D; Video V4) and EF3 with GRK1-A11 (Figure 3E; Video V5) suggested that each EF-hand independently communicates with the target, with no overlapping of the pathways.

Overall, PSN pathway analysis indicated that intermolecular information transfer concerning Ca²⁺-binding was routed by selective pathways connecting EF-hands with GRK1.

The connectivity of a PSN is an interesting feature that can be used to compare different states of the same protein and to identify hub residues (Tables S1–S3), that is amino acids involved in persistent interactions with many other residues, therefore playing a key role in maintaining network topology (Raimondi et al., 2011; Fanelli et al., 2013; Marino and Dell'orco, 2016).

An example of PSN is represented in Figure 4, where all persistent interactions in Ca²⁺-loaded NCS1 and Rec are shown. Although the number of connections (and consequently hubs) identified in NCS1 and Rec are different due to intrinsic protein properties, the highest degree hubs (degree = 8) in both proteins share their function, as NCS1 D73 and Rec D74 are the first Ca²⁺-coordinating residues of EF2 Ca²⁺-binding loop, while NCS1 D109 and Rec D110 are the equivalent for EF3.

Hub analysis was performed for all NCS1 (Table S1) and Rec (Tables S2, S3) states and was compared to our previously published data (Marino and Dell'orco, 2016) on EF3-Mg²⁺ and Ca²⁺-loaded human GCAP1 (Table S5, and Table S4 in the reference).

Multiple Sequence Alignment (Figure 5) showed that these three NCS proteins have 9 identical and 6 conserved hubs whose maximum degree was at least 5 in any of the simulated states. Interestingly, 4 identical hubs represent the first and the last residue of the Ca²⁺-binding loops of EF2 and EF3, implying that any mutation of these residues could have a detrimental effect on protein structure/function relationship, as previously demonstrated (Kitiratschky et al., 2009; Dell'orco et al., 2010; Vocke et al., 2017; Marino et al., 2018).

Finally, multiple sequence alignment was run for each protein with their respective UniRef50 cluster (proteins sharing at least 50% s.i. and 80% sequence coverage), and the conservation of each hub residue among homolog proteins was calculated (Tables S1–S3, S5) and compared to their degree (Figure 6). While no apparent mathematical relationship could describe this behavior in general, a clear, protein-specific trend could be identified. All three proteins reported a decrease of the average conservation concurrent with the decrease of the hub degree, with a 62% probability to be linearly correlated. Although the R² value of 0.62 is not sufficiently high to clearly infer a correlation between the average conservation and the hub degree, it has to be noticed that there is an intrinsic limitation in the number of points in the x-axis, as the hub degree is a discreet quantity. Therefore, R² value is greatly affected by the presence of some poorly conserved (<50%) high degree hub residues of Rec, namely I155, V87, E181, and K194.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.