Increased functional unit flexibility and solvent accessibility favours oxygen capture in molluscan hemocyanin

Abstract

Hemocyanins are a class of copper-based oxygen transport proteins, widely studied for their unique oxygen-binding processes and their role in the molluscan immune response. In this study, we utilised computational simulations to investigate the first functional unit (FU-a) of Crepidula fornicata (slipper limpet) hemocyanin, a member of the keyhole limpet hemocyanin family. Using quantum mechanics/molecular mechanics (QM/MM) methods, we designed oxygenated and deoxygenated models of FU-a and conducted molecular dynamics simulations to explore their functional dynamics and oxygen accessibility. We specifically focused on understanding the global and localised dynamics between the two conformational states. By employing principal component analysis (PCA) and modevector analysis, we differentiated the dynamic properties of the deoxygenated and oxygenated states of the hemocyanin. Furthermore, we explored the impact of oxygenation on hydration and tunnel cavity formation. Our results reveal that oxygen entry is mediated by a single bidirectional tunnel, with its permeability tightly regulated by differential histidine-based copper coordination. Importantly, we identified Glu352 as an evolutionary conserved molecular “shutter,” whose conformational changes govern the opening and closure of this tunnel. These findings provide insight into the mechanistic regulation of oxygen transport in molluscan hemocyanins, with implications for understanding their functional versatility and potential applications.

---

Introduction

Hemocyanins are large copper-containing multimeric oxygen transport proteins found in arthropods and molluscs.¹ In addition to their role in oxygen transport, hemocyanins have been implicated in the molluscan immune response, both as phenoloxidase enzymes, where proteolysis or partial denaturation enables tyrosine and derivatives to access the copper–oxygen centre,^2,3 and as a source of antimicrobial peptides.^4,5 These functional roles underline the importance of understanding the 3D-structure and dynamics of molluscan hemocyanins, including their subunit arrangement and interfaces, structural transitions and oxygen binding characteristics.

While hemocyanins from arthropods are hexamers,⁶ molluscan hemocyanins are typically either decamers or didecamers comprised of 10–20 subunits.¹ The most widely studied of the molluscan hemocyanins, the keyhole limpet class, are typically didecamers of 20 subunits arranged with d5 symmetry, where each decamer is made up of 5 dimer ‘plates’ arranged in 5-fold rotational symmetry.⁷ Each dimer plate contains 2 identical subunits but is effectively a heterodimer due to the differential structural arrangement of functional units (FUs) within the two monomers.⁸ Each monomer, approximately 3500 amino acids in length, is composed of 8 FUs (FU-a–FU-h), which contain an alpha-helical core domain containing a Type 3 copper–oxygen centre, and an adjacent beta-structure thought to regulate access to the copper–oxygen centre.⁹ A molluscan hemocyanin didecamer therefore contains 160 oxygen binding sites, which also show cooperativity,² likely arising from interactions between FUs from the same or different subunits,^10,11 although the precise molecular mechanisms remain to be elucidated.

Oxygen binding to hemocyanins has been extensively studied, and the copper–oxygen centres of each FU are highly conserved across arthropods and molluscs, with 2 copper ions coordinated to 6 histidine residues.⁷ In the deoxygenated form, the Cu ions are in the Cu(I) oxidation state, whereas in the oxygenated form (Cu₂O₂), they are oxidised to the Cu(II) state. Crystallographic studies of horseshoe crab (Limulus polyphemus) hemocyanin have demonstrated that upon oxygen binding, the Cu–Cu distance shortens by approximately 1 Å.^12,13 The oxygen binding structural transition in horseshoe crab hemocyanin has also been studied using molecular dynamics (MD) simulations, investigating the effect of hydration on facilitation of oxygen access to the copper centres,¹⁴ and how inter-subunit dynamics controls the formation of tunnels for oxygen in the hemocyanin structure.¹⁵ Most recently, a quantum mechanics/molecular mechanics approach (QM/MM) has been employed to probe the oxygenation process in detail, showing the formation of a stable superoxide intermediate consistent with previous experimental data.¹⁶

In this work we set out to analyse the molecular dynamics of oxygen binding in a molluscan hemocyanin of the keyhole limpet class, slipper limpet (Crepidula fornicata) hemocyanin (SLH), building on our recent work determining its structure using a combination of cryoEM and Alphafold2.¹⁷ We modelled the 3D-structure of FU-a, generated custom parameters for the Cu₂ and Cu₂–O₂ deoxygenated and oxygenated states, and show that oxygen entry is mediated by a single bidirectional tunnel, with its permeability tightly regulated by differential histidine-based copper coordination. We identify Glu352 as a conserved molecular “shutter,” whose conformational changes govern the opening and closure of this tunnel. We also demonstrate that the deoxygenated state displays increased flexibility and solvent accessibility, proposing that this favours oxygen capture.

Methods

Starting structures

Slipper limpet hemocyanin naturally forms large multimeric assemblies, typically decameric and didecameric in size.¹⁷ The sequence for FU-a was extracted from the SLH1 protein sequence (GenBank accession number WCA44164.1), its signal sequence was removed and the 3D structure (deoxygenated and devoid of copper ions) modelled using Alphafold3.^17,18 To accurately model the copper–oxygen centre, AlphaFill was employed to incorporate Cu₂ clusters into both the deoxygenated (FU-a_deoxy) and oxygenated (FU-a_oxy) states, with the latter additionally incorporating oxygen.¹⁹ The resulting coordination geometries were carefully analysed alongside our previously developed, unpublished AlphaFold2-based model of SLH, yielding a root-mean-square deviation (RMSD) of <1 Å. Following this, the structures were further refined using the following QM/MM methods to achieve optimal electronic and structural fidelity.

Computational methods

Molecular dynamics (MD) simulations were conducted using the GROMACS software package.²⁰ Due to the absence of force field parameters for the protein's coordination of the two metal sites, custom parameters were generated using the MCPB.py module of the AMBER suite.²¹ Bond angles, bond lengths, and partial charges were calculated through Gaussian09 software, with B3LYP/Def2-TZVP.^22,23 Dicopper clusters were assigned overall charges of +1 and +2, respectively, corresponding to the FU-a_deoxy and FU-a_oxy states, both with singlet multiplicities. Harmonic force constants and atomic charges were derived from DFT calculations using the Seminario method and restrained electrostatic potential (RESP) fitting scheme.^24–27 The resulting parameters were then integrated with the AMBER's 2019 Stony Brook forcefield (ff19SB) using the LEaP module.²⁸ Protonation states of titratable residues were determined at pH 8.07 and a salt concentration of 0.6 mM, simulating seawater conditions typical of the North Atlantic Ocean.^29,30 The systems were then solvated using the explicit water model TIP3P.³¹ Temperature coupling was applied using the v-rescale thermostat, while Particle-Mesh Ewald (PME) was employed for long-range electrostatics. Simulations were performed under the NPT ensemble with periodic boundary conditions (PBC), at a temperature of 285.43 K and a pressure of 1 atm.³² After energy minimisation and equilibration, five independent molecular dynamics simulations were executed for each system, consisting of 500 million steps with a 2 fs timestep, yielding a total simulation time of 5 microseconds per system.

Trajectory analysis

Root-mean-square deviation (RMSD), root-mean-square fluctuations (RMSF), and radius of gyration (R_g) were monitored to assess protein stability. Principal component analysis (PCA) was employed to extract dominant conformational motions and delineate structural differences between states. The covariance matrix of atomic fluctuations was constructed and diagonalised, yielding eigenvalues and corresponding eigenvectors, with eigenvalues normalised and ranked to prioritise principal motions. Eigenvectors were subsequently extracted, normalised, and sorted by their significance. A dot product matrix was computed for the leading principal components, quantifying subspace overlap and characterising shared conformational dynamics between states. Complementary metrics, including subspace overlap and maximal dot products, provided a comprehensive measure of similarity in dominant motions. Gaussian kernel density estimation was applied to principal component projections, generating contour maps that illustrate the distribution of sampled conformations, offering both quantitative and visual insights into the structural landscapes explored throughout the simulations.

Binding free energy calculation

Binding free energy was computed using gmx_MMPBSA, selected for its compatibility with GROMACS trajectories and its capacity to evaluate both Poisson–Boltzmann Surface Area (PBSA) and Generalized-Born Surface Area (GBSA) methods.³³ Calculations were performed on the equilibrated section of the MD simulation (100–1000 ns), focusing on the energetics between coordinating histidines and the active-site cargo. Given the significant electrostatic polarisation within the metal center, the internal dielectric constant was varied between 1.0 and 5.0 to assess its impact on the free energy calculations. For GBSA, a GB-Neck2 model was employed to accurately capture complex solvation environment presented in this study.³⁴ PBSA calculations followed similar parameters; however, given that the GB-Neck2 model was too computationally expensive with the available resources, the GB-OBC1 model was used instead.³⁵ In all cases, the ΔG_binding was reported as relative binding free energy, with the contribution of conformational entropy excluded from the evaluation.

Tunnel formation analysis

Dynamic tunnel formation was analysed using CAVER 3.0.³⁶ Every tenth frame from the equilibrated segment (100–1000 ns) of each replicate was extracted, resulting in a total of 90 000 frames for analysis encompassing both FU-a_deoxy and FU-a_oxy. These frames were processed following the CAVER protocol, with copper ions designated as reference points for tunnel identification. To ensure the identification of only sufficiently sized tunnels, a probe radius of 1.5 Å, corresponding to the van der Waals radius of transported oxygen, was applied. A clustering threshold of 3.0 Å was employed to refine the accuracy of tunnel detection. The bottleneck radius was then calculated to determine the maximum probe size that could pass through the narrowest section of each tunnel. Tunnels were further assessed based on throughput, ranging from 0 to 1, where a higher value indicated greater functional relevance of the channel.³⁶

Residue interaction network analysis

Residue interaction networks (RIN) were constructed to examine the alterations in local interaction networks upon oxygen binding. Analysis focused on the equilibrated segment of each replicate (100–1000 ns), with every tenth frame extracted in PDB format. The Residue Interaction Network Generator 4.0 (RING 4.0) was utilised to quantify contact frequencies across the entire protein, incorporating predefined physicochemical properties and geometric constraints in its calculation framework.³⁷ These probabilities were subsequently averaged across replicates and compared between states. RIN data was processed using in-house Python scripts, leveraging NetworkX and Matplotlib to visualise the interaction networks, where amino acids were represented as nodes, and their interactions as edges.^38,39 A key step in the analysis was identifying signal-propagating nodes within a maximum distance of 2 steps from the highlighted residue, Glu352, using a threshold of 0.3 for interaction strength. Contacts were classified as exclusive to FU-a_oxy, FU-a_deoxy, or shared between both states.

Results and discussion

Oxygen capture in molluscan hemocyanins is a highly conserved process driven by copper–oxygen interactions across multiple copper-binding sites. In deoxygenated hemocyanin, each copper ion is trigonally coordinated by three histidyl imidazole nitrogens.^{12,13,40–42} Upon oxygenation, oxygen binds as a peroxide, forming a μ:η²-η² bridge between the two copper atoms.⁴³ This binding induces a shift in oxidation state from Cu(I) to Cu(II), stabilising the peroxo-like configuration.

In addition to the peroxo form, oxygen can also bind in the bis-μ-oxo form, which is crucial for studying the flexibility and reactivity of copper–oxygen sites in hemocyanins.^44–46 In this form, the Cu–Cu distance reduces to approximately 2.8 Å, compared to 3.6 Å in the peroxo form, enhancing the interaction between copper and oxygen.⁴⁷ The isomerisation between the peroxo and bis-μ-oxo forms plays a critical role in oxygen activation and O–O bond cleavage, key processes for both oxygen binding and reactive oxygen species (ROS) deactivation in hemocyanins. Theoretical models indicate that this bis-μ-oxo form is energetically stable, with interactions between Cu 3d and O₂ π* orbitals contributing to its stability. These models suggest the bis-μ-oxo form offers a favourable pathway for oxygen species activation.^48–51 Given its stability and the pronounced difference in Cu–Cu distance, we selected the bis-μ-oxo form for our MD simulations. This choice enhances the distinction between the oxygenated and deoxygenated states, while remaining biologically relevant to hemocyanin function.

Previous studies have explored the parametrisation and subsequent simulation of horseshoe crab (arthropod) hemocyanin structures with resolved copper–oxygen sites.^14,15 We have extended this work by modelling copper–oxygen sites into slipper limpet (molluscan) hemocyanin, using a model where these sites have not yet been resolved structurally. We first aimed to optimise the orientation of the copper–oxygen site in our model with high accuracy, permitting MD simulations.

Parameterisation of metal sites

To facilitate the simulation of SLH FU-a, we first parameterised two distinct FU-a models: FU-a_oxy, containing a Cu₂O₂ core, and FU-a_deoxy, featuring a fully deoxygenated Cu₂ core (Fig. 1), using AlphaFill to incorporate the Cu₂/Cu₂O₂ site.^19,52 The resulting structures were then parameterised using the MCPB tool and subjected to geometric optimisations at the B3LYP functional level with the Def2-TZVP basis set.^21,23 Geometrical parameters, including bond lengths and angles for the FU-a active sites, are presented in Tables S1 and S2.† QM calculations reveal that in FU-a_deoxy, copper coordination forces are distributed relatively uniformly among the coordinating histidines, with minimal deviations in equilibrium bond lengths and force constants across both histidine triplets. In contrast, for FU-a_oxy, a distinct coordination pattern emerges: two histidine ligands exhibit strengthened interactions with the copper centre, while the third shows weakened association. This asymmetric coordination is replicated at both coordination sites, and results in a reduced interaction involving His71 and His215. This suggests a potential local destabilisation or altered electronic environment at these positions, which could affect overall structural stability, or influence reactivity in the Cu₂ core.

Fig. 1 Structural representation of the FU-a core highlighting the Cu₂O₂ core and imidazole groups from the coordinating histidines. Circled areas represent minimum energy structures obtained at the B3LYP/Def2-TZVP level of theory of the oxygenated (FU-a_oxy; top) and deoxygenated (FU-a_deoxy; bottom) core. Atom colour codes: Cu (orange), O (red), C (green), N (blue), H (white).

The stability of both the coordination models was then evaluated by computing the root-mean-square deviation (RMSD) of activate site heavy atoms over time. Mean RMSD values, averaged across the equilibrated region of each MD replicate, indicate greater stability in FU-a_oxy (0.52Å ± 0.11) relative to FU-a_deoxy (1.00Å ± 0.19). This suggests that oxygen binding enhances rigidity within the active site, allowing ligating residues to adopt a more stable conformation more rapidly than in the deoxygenated state.

Additionally, the minor increase in RMSD variability observed in FU-a_deoxy replicates is consistent with the increased unoccupied volume within the active dicopper centre, permitting greater conformational flexibility among ligating residues. To establish which ligating residues contributed most to this difference in coordination stability, equilibrium bond lengths obtained from QM-optimised active sites were compared to those derived from MD trajectories. We observe a mean bond length increase of 1.9% in FU-a_oxy and 12.6% in FU-a_deoxy, with His73 and His62 exhibiting the most substantial elongation for each state respectively (Table S3†). While these changes deviate from the ideal bond lengths predicted by DFT calculations, the standard deviation of the copper-ligand distances remains below 0.15 Å in all cases. This suggests that, despite the slight elongation, the bonds retain a high degree of stability, preserving the active coordination environment of the copper centers.

Moreover, when comparing the inter-copper distance (d_Cu–Cu), a distinct change occurs upon oxygen binding, reducing d_Cu–Cu from 5.07 Å in FU-a_deoxy to 2.84 Å in FU-a_oxy – consistent with the copper association observed in other dicopper–oxygen complexes.^13,53–55 Overall, these results support the stabilising influence of oxygen by way of elevated rigidity within the active site, a concept previously reported in Limulus polyphemus.¹⁴ The observed bond distances accurately reflect the conserved transformation undertaken by the dicopper centre upon oxygenation. While slight elongation of ligating distances is evident, their stability remain intact. This elongation likely reflects the influence of dynamic simulation conditions, which capture temperature, pressure, and solvent effects not represented in prior theoretical studies.

Active site energetics

Building on the structural observations discussed previously, we next evaluated the energetics of active site coordination in both oxygenated and deoxygenated states of FU-a. Specifically, we quantified the binding free energy (ΔG_bind) between the coordinating histidines and the active site cargo: Cu₂ (in the case of FU-a_deoxy) and Cu₂O₂ (for FU-a_oxy) using both MM/GBSA and MM/PBSA end-state methods (Table 1). Accurate implementation of the gmx_MMPBSA tool requires careful parameter definition, particularly the internal dielectric constant (IDC), to reliably assess binding free energies. To enhance the robustness of our approach, as such, we broadened our IDC parameter range, drawing on insights from a previous study that assessed the influence of parameter selection on calculation accuracy.⁵⁶

Table 1 Binding free energies from the MM/GBSA and MM/PBSA calculations for FU-a_oxy and FU-a_deoxy systems

ΔGbind (kcal mol^−1)	FU-aoxy	SD	FU-adeoxy	SD
MM/GBSA	−17.72	1.35	−1.29	3.04
MM/PBSA	−19.85	1.42	−4.78	3.01

---

In FU-a_deoxy systems, ΔG_bind values were largely insensitive to variations in the internal dielectric constant (IDC) (Tables S4 and S5†). In contrast, FU-a_oxy systems exhibited pronounced sensitivity, with ΔG_bind values ranging from −72.06 to −17.50 kcal mol⁻¹ using MM/GBSA and from −80.31 to −19.85 kcal mol⁻¹ using MM/PBSA. Although increasing the IDC resulted in progressively less favourable binding energies, it concurrently reduced variability across replicates. This trend suggests that incorporating greater electronic polarisation improves the stability and convergence of free energy estimates. Although the dependence of binding energy on dielectric constant used indicates that absolute binding energy from MM/PBSA approach may not be reliable, as discussed by Roux and Chipot, comparison of oxygenated with deoxygenated form should allow us to qualitatively determine relative binding energy.⁵⁷ Consequently, an IDC value of ε = 5 was selected for both states to ensure consistency and methodological robustness.

Comparison of the final free energy values highlights the notable influence of oxygenation on active site energetics. MM/GBSA calculations revealed a ΔG_bind difference of 16.43 kcal mol⁻¹ between FU-a_deoxy and FU-a_oxy, with the latter displaying more favourable binding (Table 1). Similarly, MM/PBSA calculations indicated a 15.07 kcal mol⁻¹ difference, reinforcing the enhanced thermodynamic stability of the oxygenated coordination environment. These findings are consistent with the proposed compaction and increased rigidity of the active site upon oxygen binding and underscore the intrinsic energetic favourability of oxygen coordination – an expected feature for proteins of this nature.

This energetic preference is consistent with prior studies of type-3 copper containing metalloproteins. Saito et al. reported a binding energy of −18.9 kcal mol⁻¹ for oxygenation in the T-state of Limulus polyphemus hemocyanin, a value that overestimates the experimental binding enthalpy range of −11.5 to −6.0 kcal mol⁻¹.^58–60 In contrast, a more modest binding energy of −2.3 kcal mol⁻¹ was reported by Fasulo et al., who used the NEVPT2 method to obtain a result which better aligns with experimental observation.¹⁶ The 15.07 kcal mol⁻¹ difference observed in this study falls within the expected range of hemocyanin-like systems, particularly given the use of classical MD-based free energy methods, in contrast to the QM or QM/MM approaches used in previous work, alongside expected differences in protein environment. Nonetheless, the consistent stabilisation observed upon oxygen binding in FU-a_oxy strengthens our confidence in the present parametrisation, and mirrors trends reported across both theoretical and biological systems.

Global dynamics

To gain a deeper insight into the functional consequence of oxygenation, we extended our analysis beyond the active site to examine the broader dynamical features of FU-a. To ensure an accurate representation of core dynamics, only residues 1–407 were included, thereby excluding the highly flexible C-terminal tail which would otherwise link adjacent functional units in higher order hemocyanin assemblies. RMSD and radius of gyration (R_g) were used as real-time indicators of global conformational stability and compactness (Fig. 2A and B). Mean RMSD values calculated over the equilibrated portion of the trajectories were 2.10 Å ± 0.18 for FU-a_oxy and 2.46 Å ± 0.21 for FU-a_deoxy, indicating a modest increase in global stability upon oxygen binding. This trend was reaffirmed by R_g analysis, which revealed a slight increase in compaction of FU-a_oxy (21.33 ± 0.07 Å) relative to FU-a_deoxy (21.51 ± 0.11 Å). Collectively, these observations suggest that oxygenation enhances structural compaction and promotes a more ordered global protein architecture, a phenomenon that has also been reported in previous studies of higher order hemocyanin assemblies.^14,15,61

Fig. 2 A comparison of global structural stability and residue-level flexibility between FU-a_oxy and FU-a_deoxy replicates. (A) Average Ca RMSD values for FU-a_oxy (red) and FU-a_deoxy (blue) calculated over the equilibrated portion of the trajectories. (B) Average R_g for FU-a_oxy (red) and FU-a_deoxy, indicating a modest increase in structural compaction upon oxygenation. (C) Front and rear views of FU-a_oxy, with residues exhibiting RMSF values of ≥ 1 Å highlighted in red. (D) Corresponding views of FU-a_deoxy, with highly mobile residues highlighted in blue.

Root-mean-square fluctuation (RMSF) was employed to assess residue-level flexibility across the equilibrated regions of the simulation and revealed broadly similar patterns between FU-a_oxy and FU-a_deoxy (Fig. 2C and D). To better resolve subtle localised differences, per-residue values for FU-a_oxy were subtracted from those of FU-a_deoxy. Residues exhibiting absolute differences exceeding ± 0.5 Å were subject to further analysis (Fig. S1†). As anticipated, the largest discrepancies in flexibility were predominantly localised to loop regions, specifically residues Pro46–Ala56, His133–His137, Ser153–Pro154, Ala350, and His374–Phe378. Among these, the loop spanning Pro46–Ala56 was particularly noteworthy, given its apparent interactions with neighbouring lops near the active site. The increased mobility observed in the deoxygenated state suggests a loss of inter-loop connectivity, potentially destabilising the local architecture. Interestingly, His43 and His62 – key residues involved in copper coordination – are positioned near regions of elevated flexibility and are closely spaced in the primary sequence. This suggests a potential link between the local structural dynamics and the coordination of the dicopper center. Together with the observed dynamic shifts in these loops and their proximity to the active site, these findings support a possible role for these regions in facilitating oxygen entry from the solvent environment.

Principal component analysis

While the analysis above provides valuable insight into structural stability and local fluctuation, they do not capture the collective, correlated motions that often underpin functional dynamics. To explore these broader conformational patterns and directly compare dominant motion profiles between the oxygenated and deoxygenated states, we performed principal component analysis (PCA). This dimensionality reduction technique allowed us to identify the essential dynamical subspaces by isolating the most pronounced motions sampled during the simulation.

Initial assessment of the variance distribution via scree plot (Fig. 3A) revealed that the first principal component (PC1) accounted for a disproportionately large portion of the variance in both FU-a_oxy and FU-a_deoxy. This finding suggested that PC1 primarily reflected motions of the highly flexible C-terminal tail, a region excluded from previous analysis due to its intrinsic disorder. As such, PC1 was excluded from subsequent evaluation to prevent it from dominating the interpretation of the collective dynamics. This adjustment increased the number of components required to explain 95% of the total variance – from 13 to 28 in FU-a_oxy, and from 5 to 24 in FU-a_deoxy – offering a more accurate view of the dynamics involved in the oxygenation process.

Fig. 3 Principal component analysis (PCA) of the collective motions in FU-a_oxy and FU-a_deoxy. (A) Scree plot showing the distribution of variance explained by the top 20 PCs for both oxygenated and deoxygenated states. (B) Subspace comparison based on dot product analysis of eigenvectors. (C) Projection of conformational sampling along PC2 and PC3 visualised as density heatmaps. (D) Structural visualisation of the dominant collective motions represented by PC2 in both FU-a_oxy and FU-a_deoxy. Modevectors represented as arrows depict the direction and magnitude of atomic displacement.

To further investigate the dynamic similarity between the oxygenated and deoxygenated states, a pairwise comparison was performed using the dot product matrix of eigenvectors derived from PCA. This analysis quantitively assessed the overlap in dynamical subspace between FU-a_oxy and FU-a_deoxy, revealing a subspace overlap of approximately 47.8%. This result aligns closely with the 45.7% overlap reported in Bux et al. (2021), indicating a degree of dynamic similarity between the two systems. However, some divergence from the reported value is expected due to species-specific variations in hemocyanin, as well as the order of assembly analysed.¹⁵

To directly compare the nature of the collective motions between the oxygenated and deoxygenated states, we calculated the dot product matrix of eigenvectors derived from the PCA. The analysis revealed an average maximum dot product if 0.774, indicating a substantial degree of localised similarity across several dominant PCs, suggesting notable overlap in the collective motions between the two states. The highest similarity was observed between FU-a_deoxy_PC2 and FU-a_oxy_PC2, with a dot product of 0.878, highlighting a particularly well-conserved mode of motion between the oxygenated and deoxygenated states (Fig. 3B and C).

Despite the overall similarity, subtle yet meaningful divergence was observed between the two PCs identified above to be most similar. To explore this further, we extracted the extreme conformations associated with FU-a_deoxy_PC2 and FU-a_oxy_PC2 and visualised their motions using directional modevector analysis. This method enables the simultaneous assessment of both the direction and magnitude of atomic displacement from average structures obtained from MD replicates. While most structural changes occurring in FU-a_oxy appeared conserved, though more pronounced, in the deoxygenated state, a notable exception occurred in which a loop region adjacent to the active site exhibited a high degree of displacement, a feature that was distinctly absent in FU-a_oxy replicates (Fig. 3D).

This specific loop had previously been highlighted in RMSF analysis, where it demonstrated a heightened flexibility in the deoxygenated state, an effect attributed to disruption coupling with an adjacent loop (Fig. S1†). Taken together, the findings from PCA and modevector analysis thus reinforce earlier observations, suggesting that oxygen binding plays a key role in global stabilisation, with particular impact on the loop spanning residues Pro46–Ala56 where the most pronounced conformational changes were observed.

Active site accessibility

Building on the observation that oxygen binding induces distinct collective motions, specifically within loop regions proximal to the active site, we next assessed whether these changes in conformational dynamics were accompanied by changes in solvent accessibility. Solvent accessible surface area (SASA) was computed across each of the FU-a_oxy and FU-a_deoxy ensembles to quantify differences in solvent exposure attributable to oxygenation state (Fig. 4A). Globally, the FU-a_deoxy exhibited a higher SASA (1829.96 Ų ± 42.92) relative to FU-a_oxy (1772.91 Ų ± 28.12). Although the absolute difference was modest, the temporal consistency of SASA values between groups suggests a more compact and solvent-shielded architecture in FU-a_oxy, consistent with the enhanced structural stability observed in previous analyses.

Fig. 4 Solvent accessibility analysis of FU-a_oxy and FU-a_deoxy. (A) Temporal profile of total solvent-accessible surface area (SASA) for FU-a_oxy (red) and FU-a_deoxy (blue). (B) Per-residue differential SASA calculated as SASA_deoxy – SASA_oxy. (C) Structural mapping of per-residue solvent accessibility changes, with the proposed solvent cavity explicitly labelled. Residues are coloured on a gradient scale from blue (highest solvent accessibility) to red (lowest), illustrating the spatial distribution of exposure across the structure.

To elucidate the localised effects of oxygenation on solvent accessibility, residue-specific SASA differences were calculated between the FU-a_oxy and FU-a_deoxy ensembles (Fig. 4B). Notably, Gln122 and Thr358 exhibited notable reductions in solvent exposure upon oxygen binding (ΔSASA ≤ −0.2 nm²), while Pro46, Lys119, Glu203, Asn351, and Gln376 demonstrated marked increases (ΔSASA ≥ +0.2 nm²). Structural mapping of these differences revealed a spatially coherent pattern, with pronounced increases in solvent accessibility concentrated around the active site (Fig. 4C). Among these, Pro46, located on a loop previously identified as exhibiting the greatest conformational divergence between FU-a_oxy and FU-a_deoxy, showed a substantial increase in exposure. Despite this, neighbouring residues within the same loop remained largely unaffected, indicating that the loop may not serve as a general gate for oxygen entry. Instead, these findings suggest that discrete residues within the active site microenvironment may play a more direct role in facilitating oxygen access.

These findings are consistent with prior theoretical investigations that identified solvent-accessible cavities at subunit interfaces within hemocyanin hexamers.^12–15 Notably, Bux et al. (2021) reported a pronounced increase in solvent accessibility in the Limulus polyphemus hexamer upon deoxygenation, attributing this change to the collapse of solvent tunnels following oxygen binding.¹⁵ In parallel, our data reveal localised increases in solvent exposure near the active site, reinforcing the notion that dynamic modulation of solvent-accessible tunnels is a key component of the oxygenation mechanism. The convergence of these observations underscores the importance of conformational flexibility and transient tunnel formation in regulating oxygen transport within hemocyanin.

Domain rotation

Together, the PCA and modevector analyses revealed both conserved and state-specific collective motions, underscoring the critical role of flexibility in the oxygenation mechanism. To further explore functional transitions at the quaternary level, we next focussed on domain rotation, a hallmark conformational response in hemocyanins linked to oxygen-mediated allosteric signalling. Structural comparison of Limulus polyphemus and Panulirus interruptus identifies a 7.5° rotation of domain 1 relative to domains 2 and 3, suggesting a concerted shift in the quaternary assembly upon oxygen binding.¹² Similarly, analysis of Limulus polyphemus subunit II reveals a 3.2° rotation between the two trimers within the hexamer, reflecting a rotation at the subunit interface rather than at the domain level.¹³ This rotation is thought to be central to the regulation of cooperativity in higher order assemblies, where the movement of domain 1 induces an allosteric signal which propagates throughout the entire structure, thereby enhancing oxygen affinity cross all functional units. While this transition is well-characterised in arthropod hexamers, species-specific variation may suggest that the mechanism is not universally conserved. In our model, each functional unit consists of two domains – domain I (residues 1–305), predominantly helical, and domain II (residues 306–423), characterised by its β-sheet-rich architecture – rather than the previously described 3 domain structure. These structural differences raise the question of how the allosteric mechanism may vary across species. To address this, we quantified domain rotation in our FU-a models by aligning every 10th frame of domain I to a single FU-a_deoxy reference structure. Centre-of-mass calculations were performed to determine the relative movement of domain II with respect to the fixed domain I (Fig. S2†). This analysis revealed a 1.46° inward rotation of domain II towards the active site upon oxygen binding, a motion that was deemed to be highly statistically significant (P < 0.0001) (Fig. S3A†).

Intriguingly, Magnus et al. and Hazes et al. both suggest that domain rotation directly impacts the interaction between Phe49 and coordinating histidines within the dinuclear copper site. Specifically, they describe how the aromatic ring of Phe49 packs tightly against the imidazole ring of His328 in the closed state, restricting His328's movement and preventing it from adopting the optimal position required for oxygen coordination. This restriction is lifted in the open conformation, allowing His328 to adopt a lower energy position, thereby facilitating increased oxygen affinity.^12,13 However, our network analysis reveals a divergent behaviour in slipper limpet hemocyanin. Specifically, in the FU-a_oxy state, the corresponding residue, Phe67, forms persistent interactions with coordinating histidines His62, His184, and His215. Notably, the incidence of these interactions decreases by 40.3%, 46.1%, and 47.0%, respectively, upon deoxygenation (Fig. S3B†). This suggests that, in contrast to other species, Phe67 does not obstruct the histidines from adopting the energetically favourable conformation required for optimal oxygen binding. Rather it appears to play a role in stabilising the dicopper site in the oxygenated state, supporting the structural integrity of the active site without hindering oxygen coordination.

Tunnel formation

To better understand how solvent accessibility near the active site varies between functional states, we conducted tunnel detection analysis to identify transient pathways that could facilitate oxygen diffusion. A summary of the detected tunnels is presented in Table 2, with the top-ranked cluster identified as the most probable candidate based on its higher throughput and average bottleneck radius. The analysis revealed a pronounced disparity in tunnel formation between the two functional states. The FU-a_deoxy system exhibited a 13.6-fold increase in tunnel frequency, with 776 tunnels detected across the simulation ensemble, compared to just 57 tunnels in FU-a_oxy. This highlights a clear bias in tunnel permeability towards the deoxygenated state.

Table 2 Tunnel analysis of FU-a_oxy and FU-a_deoxy states performed using Caver 3.0 with a probe radius of 1.5 Å and clustering threshold of 3.0 Å. Bottleneck radius (BR) and throughput (TP) values describe tunnel size and accessibility

ID	No_snaps	Avg_BR	SD	Max_BR	Avg_L	SD	Avg_thru	SD
FU-aoxy
1	27.80	1.53	0.03	1.62	15.94	1.50	0.63	0.03
2	24.20	1.53	0.03	1.59	20.57	1.61	0.56	0.03
3	17.00	1.55	0.04	1.63	24.79	2.84	0.51	0.05
4	14.25	1.53	0.02	1.57	29.69	4.83	0.46	0.04
5	7.00	1.55	0.01	1.58	37.01	1.48	0.40	0.01
FU-adeoxy
1	980.00	1.58	0.07	1.90	17.24	2.79	0.61	0.05
2	875.75	1.57	0.06	1.85	17.90	2.29	0.60	0.04
3	696.25	1.56	0.06	1.82	23.31	2.80	0.53	0.04
4	599.75	1.57	0.06	1.85	26.02	3.11	0.53	0.05
5	458.75	1.55	0.05	1.77	34.88	3.73	0.45	0.05

---

Bottleneck analysis, which quantifies the narrowest point along the tunnel, showed comparable mean values for each functional states: 1.54 Å for FU-a_oxy and 1.57 Å for FU-a_deoxy, both of which are theoretically sufficient for oxygen transport. However, in FU-a_deoxy, the bottleneck radii were observed to be more variable, with some tunnels exhibiting radii exceeding 1.9 Å. In contrast FU-a_oxy showed little deviation from the mean, with the maximum bottleneck radius observed at 1.62 Å. These results, characterised by a notable increase in tunnel formation frequency and larger cavities in the deoxygenated state, underscores the inherent favourability of FU-a_deoxy for oxygen import relative to FU-a_oxy (Fig. 5).

Fig. 5 Tunnel analysis in FU-a_oxy and FU-a_deoxy. The highest-throughput tunnel clusters are identified and shown as red spheres in (A) FU-a_oxy and blue spheres in (B) FU-a_deoxy. Minimum bottleneck radius was defined at 1.5 Å.

Given the crucial role of channel bottlenecks in controlling access to the active site, we next focused on identifying the residues that contribute to bottleneck formation and their variability across the functional states. By analysing the top-ranked tunnel clusters for both FU-a_oxy and FU-a_deoxy, we were able to pinpoint residues that consistently participate in bottleneck constriction, highlighting those most influential in regulating tunnel permeability. Residues that appeared in ≥ 70% of tunnel snapshots were considered crucial, providing a robust basis for identifying key structural determinants of channel accessibility.

The analysis revealed distinct bottleneck profiles between the two functional states. In FU-a_oxy, where the tunnel adopts a more constricted conformation, the bottleneck region is primarily shaped by residues His62, Leu202, Leu347, and Glu352. Contrastingly, FU-a_deoxy displayed a substantially altered residue composition at the bottleneck. While His62 and Leu202 remained involved in both states, the contributions of Leu347 and Glu352 were markedly reduced in the deoxygenated state. Notably, Trp355 emerged as a key residue for tunnel constriction in FU-a_deoxy, contributing to 71.11% of observed tunnel snapshots.

The consistent involvement of His62 and Leu202 across both states is unsurprising, given His62's role in copper coordination and Leu202's proximity to the active site. These residues, however, likely contribute to constriction in an accessory role, supporting but not directly controlling oxygen gating. In contrast, while Leu347 is positioned near the active site, its internal orientation directs the flow within the cavity rather than obstructing it, making its role in directly mediating oxygen import unlikely. Trp355, on the other hand, occupies a position where its steric bulk could block oxygen entry, with the indole sidechain large enough to potentially obstruct the channel. However, the lack of contribution in FU-a_oxy bottlenecks suggests that Trp355 is not crucial for blocking oxygen entry when oxygen is already bound, supporting the idea that its involvement is primarily linked to the deoxygenated state.

Glu352 was the final residue implicated in bottleneck formation in FU-a_oxy. Like Trp355, it resides on a dynamic loop near the active site but is located on the exterior of the protein, where its carboxyl sidechain may sterically obstruct the channel, potentially limiting oxygen entry. The contribution of Glu352 to the bottleneck is notable, highlighting a minor decrease upon deoxygenation. This reduction is consistent with the expected channel expansion in the deoxygenated state, allowing for greater variability in the residues that contribute to the bottleneck as the tunnel widens. Hazes et al. previously proposed that glutamic acid residues function as a molecular “shutter”, regulating oxygen entry in Limulus polyphemus, with particular emphasis on the critical role of the negative charge.¹² To assess the conservation of this mechanism, we performed sequence alignments of the functional units from slipper limpet hemocyanin, keyhole limpet hemocyanin, and octopus hemocyanin (FU-d only) (Table S7†). Additionally, we conducted a broader alignment encompassing all molluscan hemocyanin sequences (Table S8†). Our results indicate that Glu352 was fully conserved across all sequences (Table S6†), underscoring its critical, evolutionary conserved, role in regulating oxygen transport.

Mapping the detected tunnels onto the didecameric hemocyanin structure further contextualises their functional significance. The tunnel is positioned on the protein's solvent-exposed exterior, placing it in direct communication with the surrounding environment. Such an arrangement would optimise access for oxygen uptake and release, supporting the proposed dynamic gating mechanism inferred from residue contributions. Nevertheless, whether this tunnel trajectory and gating mechanism are conserved across all functional units within the oligomer remains unclear and warrants further investigation. Taken together, these findings establish a mechanistic link between active site dynamics, tunnel formation, and residue-specific gating in molluscan hemocyanins, providing new insights into how oxygen access is precisely regulated in response to functional state.

Driving forces of oxygen release

Our structural analyses of FU-a reveal distinct conformational differences between oxygenated and deoxygenated states, yet the exact mechanism of oxygen release remains unclear, especially given the stability of the oxygen-bound form. Existing evidence supports our findings, indicating that oxygen release in molluscan hemocyanins involves an increase in the Cu–Cu distance at the binuclear copper center.^12,13 However, these local changes are modulated by allosteric regulation within the multimeric assembly, governed by the surrounding ionic environment.^62–65

Since oxygen binding restricts flexibility and solvent accessibility within FU-a, dissociation is ultimately driven by external physiological cues — namely, chloride-mediated inhibition of domain rearrangement, acidification, and calcium-dependent stabilisation — that trigger quaternary level rearrangements and transitions.^63,66,67 These higher-order conformational dynamics enable the cooperative exchange between oxygenated and deoxygenated states. Consequently, a comprehensive understanding of oxygen release necessitates investigation at the multimeric level, as single-unit analyses cannot fully capture the mechanisms underlying allosteric regulation.

Network analysis

Given the strong evolutionary conservation and functional importance of Glu352 in regulating tunnel access, we next sought to investigate the structural pathways through which oxygen binding could modulate its conformation. To achieve this, we constructed residue interaction networks (RINs) from molecular dynamics simulations, aiming to map potential signal propagation routes from the oxygen-binding site to Glu352 and identify residues that mediate these transitions (Fig. 6).

Fig. 6 Residue interaction network (RIN) depicting short-range connectivity between coordinating histidines and Glu352 in FU-a_oxy and FU-a_deoxy. Interaction types are annotated, with the colour and style of connecting line indicating if the interaction is unique to FU-a_oxy (orange), FU-a_deoxy (blue), or shared (green). Interaction types include van der Waals (VDW), hydrogen bonds (HBONDS), pi–pi stacking (π–π) stacking, ionic interactions (IONIC), and π–cation interactions (PICATION).

RINs were generated using RING 4.0, with every 10th frame sampled from each trajectory to capture the dynamic ensemble. Interaction probabilities were then averaged across replicates to produce consensus networks for each functional state. Given the close spatial proximity of Glu352 to the active site, we hypothesised that signal transmission would likely occur via direct or near–direct interactions. Accordingly, we applied a stringent interaction strength threshold (≥0.3) and restricted pathway lengths to a maximum of two residues from Glu352. This focused strategy enabled us to isolate the most plausible short-range communication routes coupling oxygen binding to tunnel gating.

Subsequent analysis revealed distinct signal propagation networks associated with the oxygenated and deoxygenated states (Fig. 6). In FU-a_oxy, Glu352 is stably coordinated within a network dominated by active site residues His62 and Leu202, as well as loop residues Gly318 and Gly349. This core interaction network is further reinforced by secondary contacts, positioning Glu352 as a central hub within the structural framework. Upon deoxygenation, substantial rearrangements are observed: interactions linking Leu202, His62, His43, and Glu352 are disrupted, resulting in increased flexibility of Glu352, as indicated by elevated RMSF values (Fig. S1†).

To assess how Glu352 dynamics influence tunnel accessibility, we calculated the heavy-atom RMSD of its side chain relative to the starting conformation across structural ensembles of the oxygenated and deoxygenated states, and analysed rotamer distributions (Fig. S2†). In FU-a_oxy, Glu352 adopts two primary rotamers that retain a consistent orientation relative to His62, with only minor shifts in the carboxyl group (Fig. S4†). As a result, the Glu352 side chain remains conformationally constrained, and the tunnel remains blocked in both cases. In contrast, the deoxygenated state shows greater conformational variability. While the high-RMSD conformation preserves contacts with His62 and occludes the tunnel, the low-RMSD conformation shows a marked reorientation—Glu352 disengages from His62 and redirects its carboxyl group into the active site cavity, permitting tunnel opening. This conformational flexibility is absent in FU-a_oxy, indicating that oxygen binding triggers Cu–Cu association, which in turn drives the repositioning of His62. This movement of His62 is fundamental for stabilising Glu352 in a conformation that blocks the tunnel (Fig. S4†). Deoxygenation disrupts this coordination network, allowing Glu352 to adopt alternative conformations that transiently open the tunnel, facilitating oxygen diffusion into the active site.

To assess the conservation of this mechanism, we revisited sequence alignments performed across slipper limpet hemocyanin, keyhole limpet hemocyanin, and octopus hemocyanin (FU-d only), as well as a broader dataset of molluscan hemocyanins (Tables S7 and S8†). Unsurprisingly, key components of the rigid coordination network including His43, His62, and Leu202, are fully conserved across all functional units and molluscan sequences, supporting their essential structural roles (Table S6†). Gly349 and Trp355 are likewise fully conserved at the broader molluscan level; however, among the functional units analysed, occasional substitutions of Trp355 with phenylalanine and Gly349 with aspartate are observed. Given their respective BLOSUM62 log-odds scores of +1 and −1, these substitutions are considered conservative and unlikely to compromise network integrity. In contrast, Gly318 exhibits considerable sequence variability across both molluscan hemocyanins and functional units. This divergence suggests that, while Gly318 participates in the network, its contribution likely involves backbone-mediated rather than sidechain-specific interactions. Nevertheless, the overall high degree of conservation within this short signal propagation pathway strongly supports its functional relevance, and suggests that the mechanism of oxygen gating through Glu352-centered network rearrangement is a broadly conserved feature across molluscan hemocyanins.

Conclusions

This study presents a comprehensive investigation into the structural dynamics underpinning oxygen binding and transport in FU-a of slipper limpet hemocyanin. Our analysis revealed that oxygenation led to an overall increase in structural rigidity, a decrease in solvent accessibility and a constriction of the oxygen tunnel, suggesting that oxygen binding increases protein compactness and restricts flexibility. We identified that deoxygenation triggers conformational changes leading to an increase in the flexibility of the conserved gatekeeper residue Glu352, enhancing solvent accessibility to the active site and facilitating oxygen capture. Together, our findings provide new mechanistic insights into how hemocyanins regulate oxygen binding through dynamic modulation of tunnel formation and gating, and underscore the critical role of conserved residues such as Glu352 in the facilitation of oxygen capture.

Data availability

Data supporting this article, including pdb, topology and trajectory files for both oxy and deoxy hemocyanin are available at https://doi.org/10.5281/zenodo.15350967.

Author contributions

JGD was involved in investigation, formal analysis, methodology, writing – original draft and visualisation. JAP was involved in methodology, supervision and writing – review and editing. MTY was involved in conceptualisation, visualization and writing – original draft and review and editing. GEM was involved in conceptualisation, methodology, supervision and writing – review and editing.

Conflicts of interest

M. T. Y. has previously been supported by funding from Innovate UK (Technology Strategy Board grants 69001 and 36295). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Acknowledgements

This research was undertaken both using the supercomputing facilities at Cardiff University operated by Advanced Research Computing at Cardiff (ARCCA) on behalf of the Cardiff Supercomputing Facility and the HPC Wales and Supercomputing Wales (SCW) projects. We acknowledge the support of the latter, which is part-funded by the European Regional Development Fund (ERDF) via the Welsh Government and ARCHER granted via the UK High-End Computing Consortium for Biomolecular Simulation (http://hecbiosim.ac.uk/).

References

1. S. Kato, T. Matsui, C. Gatsogiannis and Y. Tanaka, Molluscan hemocyanin: structure, evolution, and physiology, Biophys. Rev., 2018, 10(2), 191–202.
2. H. Decker, N. Hellmann, E. Jaenicke, B. Lieb, U. Meissner and J. Markl, Minireview: Recent progress in hemocyanin research, Integr. Comp. Biol., 2007, 47(4), 631–644.
3. T. Masuda, S. Baba, K. Matsuo, S. Ito and B. Mikami, The high-resolution crystal structure of lobster hemocyanin shows its enzymatic capability as a phenoloxidase, Arch. Biochem. Biophys., 2020, 688, 108370.
4. C. J. Coates and J. Nairn, Diverse immune functions of hemocyanins, Dev. Comp. Immunol., 2014, 45(1), 43–55.
5. J. Zhuang, C. J. Coates, H. Zhu, P. Zhu, Z. Wu and L. Xie, Identification of candidate antimicrobial peptides derived from abalone hemocyanin, Dev. Comp. Immunol., 2015, 49(1), 96–102.
6. J. Markl and H. Decker, Molecular Structure of the Arthropod Hemocyanins, in, Blood and Tissue Oxygen Carriers, ed. Mangum C. P., Springer Berlin Heidelberg, Berlin, Heidelberg, 1992. pp. 325–376.
7. J. Markl, Evolution of molluscan hemocyanin structures, Biochim. Biophys. Acta, 2013, 1834(9), 1840–1852.
8. C. Gatsogiannis and J. Markl, Keyhole limpet hemocyanin: 9-A CryoEM structure and molecular model of the KLH1 didecamer reveal the interfaces and intricate topology of the 160 functional units, J. Mol. Biol., 2009, 385(3), 963–983.
9. K. N. Naresh, A. Sreekumar and S. S. Rajan, Structural insights into the interaction between molluscan hemocyanins and phenolic substrates: An in silico study using docking and molecular dynamics, J. Mol. Graph. Model., 2015, 61, 272–280.
10. M. Perbandt, E. W. Guthohrlein, W. Rypniewski, K. Idakieva, S. Stoeva and W. Voelter, et al., The structure of a functional unit from the wall of a gastropod hemocyanin offers a possible mechanism for cooperativity, Biochemistry, 2003, 42(21), 6341–6346.
11. Q. Zhang, X. Dai, Y. Cong, J. Zhang, D. H. Chen and M. T. Dougherty, et al., Cryo-EM structure of a molluscan hemocyanin suggests its allosteric mechanism, Structure, 2013, 21(4), 604–613.
12. B. Hazes, K. A. Magnus, C. Bonaventura, J. Bonaventura, Z. Dauter and K. H. Kalk, et al., Crystal structure of deoxygenated Limulus polyphemus subunit II hemocyanin at 2.18 A resolution: clues for a mechanism for allosteric regulation, Protein Sci., 1993, 2(4), 597–619.
13. K. A. Magnus, B. Hazes, H. Ton-That, C. Bonaventura, J. Bonaventura and W. G. Hol, Crystallographic analysis of oxygenated and deoxygenated states of arthropod hemocyanin shows unusual differences, Proteins, 1994, 19(4), 302–309.
14. K. Bux, S. A. Ali and S. T. Moin, Hydration facilitates oxygenation of hemocyanin: perspectives from molecular dynamics simulations, Eur. Biophys. J., 2018, 47(8), 925–938.
15. K. Bux, X. Shen, M. Tariq, J. Yin, S. T. Moin and D. Bhowmik, et al., Inter-Subunit Dynamics Controls Tunnel Formation During the Oxygenation Process in Hemocyanin Hexamers, Front. Mol. Biosci., 2021, 8, 710623.
16. F. Fasulo, A. Terán, M. Pavone and A. B. Muñoz-García, Revisiting Oxygen Transport Features of Hemocyanin with NEVPT2 Level QM/MM Calculations, J. Chem. Theory Comput., 2025, 21(4), 2108–2117.
17. G. Pasqualetto, A. Mack, E. Lewis, R. Cooper, A. Holland and U. Borucu, et al., CryoEM structure and Alphafold molecular modelling of a novel molluscan hemocyanin, PLoS One, 2023, 18(6), e0287294.
18. J. Abramson, J. Adler, J. Dunger, R. Evans, T. Green and A. Pritzel, et al., Accurate structure prediction of biomolecular interactions with AlphaFold 3, Nature, 2024, 630(8016), 493–500.
19. M. L. Hekkelman, I. De Vries, R. P. Joosten and A. Perrakis, AlphaFill: enriching AlphaFold models with ligands and cofactors, Nat. Methods, 2023, 20(2), 205–213.
20. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith and B. Hess, et al., GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX, 2015, 2, 19–25.
21. P. Li and K. M. Merz Jr, MCPB.py: A Python Based Metal Center Parameter Builder, J. Chem. Inf. Model., 2016, 56(4), 599–604.
22. M. J. Frisch, G. Trucks, H. Schlegel, G. Scuseria, M. Robb, J. Cheeseman, et al., Gaussian 16, Revision B.01, Gaussian Inc, Wallingford CT, 2009.
23. F. Weigend and R. Ahlrichs, Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy, Phys. Chem. Chem. Phys., 2005, 7(18), 3297–3305.
24. J. M. Seminario, Calculation of intramolecular force fields from second-derivative tensors, Int. J. Quant. Chem., 1996, 60(7), 1271–1277.
25. C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model, J. Phys. Chem., 1993, 97(40), 10269–10280.
26. B. H. Besler, K. M. Merz Jr and P. A. Kollman, Atomic charges derived from semiempirical methods, J. Comput. Chem., 1990, 11(4), 431–439.
27. P. Cieplak, W. D. Cornell, C. Bayly and P. A. Kollman, Application of the multimolecule and multiconformational RESP methodology to biopolymers: Charge derivation for DNA, RNA, and proteins, J. Comput. Chem., 1995, 16(11), 1357–1377.
28. C. Tian, K. Kasavajhala, K. A. A. Belfon, L. Raguette, H. Huang and A. N. Migues, et al., ff19SB: Amino-Acid-Specific Protein Backbone Parameters Trained against Quantum Mechanics Energy Surfaces in Solution, J. Chem. Theory Comput., 2020, 16(1), 528–552.
29. L.-Q. Jiang, B. R. Carter, R. A. Feely, S. K. Lauvset and A. Olsen, Surface ocean pH and buffer capacity: past, present and future, Sci. Rep., 2019, 9(1), 18624.
30. M. Zweng, J. R. Reagan, D. Seidov, T. Boyer, M. Locarnini, H. Garcia, et al., World Ocean Atlas 2018, ed. A. Salinity, Mishonov Technical, vol. 2, 2019.
31. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys., 1983, 79(2), 926–935.
32. Sea Temperature Info, https://seatemperature.info/pembroke-dock-water-temperature.html.
33. M. S. Valdes-Tresanco, M. E. Valdes-Tresanco, P. A. Valiente and E. Moreno, gmx_MMPBSA: A New Tool to Perform End-State Free Energy Calculations with GROMACS, J. Chem. Theory Comput., 2021, 17(10), 6281–6291.
34. H. Nguyen, A. Perez, S. Bermeo and C. Simmerling, Refinement of Generalized Born Implicit Solvation Parameters for Nucleic Acids and Their Complexes with Proteins, J. Chem. Theory Comput., 2015, 11(8), 3714–3728.
35. A. Onufriev, D. Bashford and D. A. Case, Exploring protein native states and large-scale conformational changes with a modified generalized born model, Proteins, 2004, 55(2), 383–394.
36. E. Chovancova, A. Pavelka, P. Benes, O. Strnad, J. Brezovsky and B. Kozlikova, et al., CAVER 3.0: a tool for the analysis of transport pathways in dynamic protein structures, PLoS Comput. Biol., 2012, 8(10), e1002708.
37. A. Del Conte, G. F. Camagni, D. Clementel, G. Minervini, A. M. Monzon and C. Ferrari, et al., RING 4.0: faster residue interaction networks with novel interaction types across over 35,000 different chemical structures, Nucleic Acids Res., 2024, 52(W1), W306–W312.
38. A. Hagberg, P. J. Swart and D. A. Schult. Exploring Network Structure, Dynamics, and Function Using NetworkX, Los Alamos National Laboratory (LANL), Los Alamos, NM, United States, 2008.
39. J. D. Hunter, Matplotlib: A 2D graphics environment, Comput. Sci. Eng., 2007, 9(03), 90–95.
40. P. Fariselli, A. Bottoni, F. Bernardi and R. Casadio, Quantum mechanical analysis of oxygenated and deoxygenated states of hemocyanin: theoretical clues for a plausible allosteric model of oxygen binding, Protein Sci., 1999, 8(7), 1546–1550.
41. K. Magnus and H. Ton-That, The crystal structure of the oxygenated form of subunit II of Limulus polyphemus hemocyanin, J. Inorg. Biochem., 1991, 43(3), 167.
42. B. Hazes and W. G. Hol, Comparison of the hemocyanin β-barrel with other greek key β-barrels: Possible importance of the “β-zipper” in protein structure and folding, Proteins: Struct., Funct., Bioinf., 1992, 12(3), 278–298.
43. S. Hirota, T. Kawahara, M. Beltramini, P. Di Muro, R. S. Magliozzo and J. Peisach, et al., Molecular basis of the Bohr effect in arthropod hemocyanin, J. Biol. Chem., 2008, 283(46), 31941–31948.
44. J. A. Halfen, S. Mahapatra, E. C. Wilkinson, S. Kaderli, V. G. Young Jr and L. Que Jr, et al., Reversible cleavage and formation of the dioxygen O-O bond within a dicopper complex, Science, 1996, 271(5254), 1397–1400.
45. J. Que and T. W. B. Lawrence, Bis(μ-oxo)dimetal “Diamond” Cores in Copper and Iron Complexes Relevant to Biocatalysis, Angew. Chem., Int. Ed., 2002, 41(7), 1114–1137.
46. S. Mahapatra, J. A. Halfen, E. C. Wilkinson, G. Pan, X. Wang and V. G. Young, et al., Structural, Spectroscopic, and Theoretical Characterization of Bis(μ-oxo)dicopper Complexes, Novel Intermediates in Copper-Mediated Dioxygen Activation, J. Am. Chem. Soc., 1996, 118(46), 11555–11574.
47. A. Bérces, Ligand Effects in the Models and Mimics of Oxyhemocyanin and Oxytyrosinase. A Density Functional Study of Reversible Dioxygen Binding and Reversible O–O Bond Cleavage, Inorg. Chem., 1997, 36(21), 4831–4837.
48. C. J. Cramer, B. A. Smith and W. B. Tolman, Ab Initio Characterization of the Isomerism between the μ-η2:η2-Peroxo- and Bis(μ-oxo)dicopper Cores, J. Am. Chem. Soc., 1996, 118(45), 11283–11287.
49. G. L. Estiú and M. C. Zerner, Calculations on the Electronic Structure and UV-visible Spectrum of Oxyhemocyanin, J. Am. Chem. Soc., 1999, 121(9), 1893–1901.
50. M. J. Henson, P. Mukherjee, D. E. Root, T. D. P. Stack and E. I. Solomon, Spectroscopic and Electronic Structural Studies of the Cu(III)2 Bis-μ-oxo Core and Its Relation to the Side-On Peroxo-Bridged Dimer, J. Am. Chem. Soc., 1999, 121(44), 10332–10345.
51. A. Berces, T. Ziegler and L. Fan, Density functional study of the harmonic force fields of cyclopentadienyl anion, cyclopentadienyl lithium, and ferrocene, J. Phys. Chem., 1994, 98(6), 1584–1595.
52. J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov and O. Ronneberger, et al., Highly accurate protein structure prediction with AlphaFold, Nature, 2021, 596(7873), 583–589.
53. M. E. Cuff, K. I. Miller, K. E. Van Holde and W. A. Hendrickson, Crystal Structure of a Functional Unit from Octopus Hemocyanin, J. Mol. Biol., 1998, 278(4), 855–870.
54. Y. Matoba, T. Kumagai, A. Yamamoto, H. Yoshitsu and M. Sugiyama, Crystallographic evidence that the dinuclear copper center of tyrosinase is flexible during catalysis, J. Biol. Chem., 2006, 281(13), 8981–8990.
55. M. Beltramini, L. Bubacco, L. Casella, G. Alzuet, M. Gullotti and B. Salvato, The oxidation of hemocyanin. Kinetics, reaction mechanism and characterization of met-hemocyanin product, Eur. J. Biochem., 1995, 232(1), 98–105.
56. M. E. Valdes-Tresanco, M. S. Valdes-Tresanco, E. Moreno and P. A. Valiente, Assessment of Different Parameters on the Accuracy of Computational Alanine Scanning of Protein–Protein Complexes with the Molecular Mechanics/Generalized Born Surface Area Method, J. Phys. Chem. B, 2023, 127(4), 944–954.
57. B. Roux and C. Chipot, Editorial Guidelines for Computational Studies of Ligand Binding Using MM/PBSA and MM/GBSA Approximations Wisely, J. Phys. Chem. B, 2024, 128(49), 12027–12029.
58. T. Saito and W. Thiel, Quantum Mechanics/Molecular Mechanics Study of Oxygen Binding in Hemocyanin, J. Phys. Chem. B, 2014, 118(19), 5034–5043.
59. A. Klarman and E. Daniel, Oxygen binding properties of stripped (calcium ion and magnesium ion free) hemocyanin from the scorpion Leirus quinquestriatus, Biochemistry, 1980, 19(23), 5176–5180.
60. K. D. Karlin, W. B. Tolman, S. Kaderli and A. D. Zuberbühler, Kinetic and thermodynamic parameters of copper-dioxygen interaction with different oxygen binding modes, J. Mol. Catal. A: Chem., 1997, 117(1), 215–222.
61. R. Sterner, T. Vogl, H. J. Hinz, F. Penz, R. Hoff and R. Foll, et al., Extreme thermostability of tarantula hemocyanin, FEBS Lett., 1995, 364(1), 9–12.
62. M. Brouwer, C. Bonaventura and J. Bonaventura, Oxygen binding by Limulus polyphemus hemocyanin: allosteric modulation by chloride ions, Biochemistry, 1977, 16(17), 3897–3902.
63. M. Brouwer, C. Bonaventura and J. Bonaventura, Chloride and pH dependence of cooperative interactions in Limulus polyphemus hemocyanin, Prog. Clin. Biol. Res., 1982, 81, 231–256.
64. M. Brouwer, C. Bonaventura and J. Bonaventura, Metal ion interactions with Limulus polyphemus and Callinectes sapidus hemocyanins: stoichiometry and structural and functional consequences of calcium(II), cadmium(II), zinc(II), and mercury(II) binding, Biochemistry, 1983, 22(20), 4713–4723.
65. J. L. Larimer and A. F. Riggs, Properties of Hemocyanins. I. The Effect of Calcium Ions on the Oxygen Equilibrium of Crayfish Hemocyanin, Comp. Biochem. Physiol., 1964, 13, 35–46.
66. B. Hazes, K. A. Magnus, K. H. Kalk, C. Bonaventura and W. G. Hol, Nitrate binding to Limulus polyphemus subunit type II hemocyanin and its functional implications, J. Mol. Biol., 1996, 262(4), 532–541.
67. A. Klarman and E. Daniel, Built-in co-operativity of oxygen binding by arthropod hemocyanins, J. Mol. Biol., 1977, 115(2), 257–261.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5ra03248b

‡ These authors contributed equally.

---