Saxitoxin and tetrodotoxin bioavailability increases in future oceans

Increasing atmospheric CO2 levels are largely absorbed by the ocean, decreasing surface water pH1. In combination with increasing ocean temperatures, these changes have been identified as a major sustainability threat to future marine life2. Interactions between marine organisms are known to depend on biomolecules, although the influence of oceanic pH on their bioavailability and functionality remains unexplored. Here we show that global change substantially impacts two ecological keystone molecules3 in the ocean, the paralytic neurotoxins saxitoxin and tetrodotoxin. Increasing temperatures and declining pH increase the abundance of their toxic forms in the water. Our geospatial global model predicts where this increased toxicity could intensify the devastating impact of harmful algal blooms, for example through an increased incidence of paralytic shellfish poisoning. Calculations of future saxitoxin toxicity levels in Alaskan clams, Saxidomus gigantea, show critical exceedance of limits safe for consumption. Our findings for saxitoxin and tetrodotoxin exemplify potential consequences of changing pH and temperature on chemicals dissolved in the sea. This reveals major implications not only for ecotoxicology, but also for chemical signals that mediate species interactions such as foraging, reproduction or predation in the ocean, with unexplored consequences for ecosystem stability and ecosystem services. Ocean warming and acidification will affect the structure and bioavailability of biomolecules. The toxic form of two neurotoxins will increase with climate change, presenting an ecotoxicology risk with global hotspots as exemplified by saxitoxin toxicity in Alaskan butter clam.

Climate change is not only increasing oceanic water temperatures, but also decreasing seawater pH as increasing atmospheric carbon dioxide (CO2) is absorbed by the ocean 1,4 . The change occurs through the formation of carbonic acid, which further dissociates into HCO3and protons, leading to a predicted drop of up to 0.4 pH units by 2100, reaching mean pH levels of 7.7 in a high-emission scenario (Representative Concentration Pathway, RCP8.5) 1 . In some coastal areas, seawater conditions of pH 7.2 are already observed temporarily 5 and are predicted to decrease further in the future. Environmental change presents a significant challenge to marine organisms at physiological, ecological, as well as behavioural level. The current rate of pH change already impacts marine organisms' calcification 6 , physiology and fitness 7 .
Interference with acid-base balance and the control of neurotransmitter function have been proposed as possible mechanisms by which ocean acidification could disrupt olfactory-mediated behaviours 8 . But there is also increasing evidence that a direct impact of pH on information-carrying signalling cues and their corresponding receptors could cause info-disruption in marine chemical communication 9,10 .
Marine organisms use a wide range of biomolecules to locate food and mating partners or to deter predators 11 . Many of these molecules also possess functional chemical groups that are sensitive to pH including hydroxyl-, carbonyl-, carboxyl-, amine-, phosphate-or sulfide-groups. Thus, changes in pH in future oceans can potentially alter a range of biological functions 2,9 . Among these, saxitoxin and tetrodotoxin have a large variety of ecological functions at very low effective concentrations 12,13 . They serve as antipredator defence through accumulation in cells, skin, tissue, eggs and oocytes in dinoflagellates, snails, ribbon worms, blue-ring octopus and pufferfish, or can be used as offensive weapon upon prey organisms 13,14 . Both toxins are released into the environment for communication purposes, e.g., as an attractant pheromone for male pufferfish 14 , or as sex-pheromone during the gametogenesis of Alexandrium sp. 15 , a dinoflagellate causing harmful algal blooms (HABs). High levels of STX usually contained within the algal cells are further released upon cell lysis at bloom-termination 16 . Predicted changes in climate are expected to further increase the duration, distribution and severity of HABs 17 , while ocean acidification has been directly shown to give toxic microalgae an advantage during a normal plankton bloom, resulting in their mass development and formation of HABs 18 . STX produced by Alexandrium sp. often accumulates in the food chain, causing paralytic shellfish poisoning (PSP) and major die-offs of fish, benthic invertebrates, and marine mammals 19,20 , with implications for marine ecosystems as well as global food security. HABs can also cause harm to humans, including direct human mortalities, mainly due to ingesting toxic seafood, direct skin contact with contaminated water or inhaling aerosolized biotoxins 21 . By blocking ion channels in nerves and muscles, STX and TTX cause partial or fatal paralysis 21 .
Both neurotoxins, STX and TTX, contain functional chemical groups that are impacted by pH 22,23 . Their protonated forms, which are more prevalent in acidified conditions, are known to possess a more effective inhibitory capacity for ion channels 24 . Once protonated, there are strong electrostatic interactions of the toxins' hydroxyl and positively charged 7,8,9-guanidinium groups with the negatively charged carboxylic side chains of the ion channels' extracellular selectivity filter site 24,25 . These lead to a full blockage of voltage-gated sodium (NaV) channels in nerves and muscles 13 , L-type Ca 2+ channels 24 and voltage-gated potassium channels 13 . In comparison to their non-or partly protonated counterparts, the fully protonated forms of STX and TTX foster an even stronger electrostatic interaction with the channels, preventing ion flux, and are therefore more potent in their toxicity 24 . The longer TTX/STX is bound, the more damaging the effect on nerve and muscle fibres 24 . The effects of pH on STX and TTX toxicity have been shown in the laboratory, but have not been translated into an ecological context, nor quantified for global future ocean models.
Here we calculate environmentally-mediated differences in protonation states of both TTX and STX and some of their derivatives in the context of published oceanic climate change scenarios 1 , visualise their toxicity-enhancing electrostatic differences and map their global abundance today and in future oceans.
We calculated the relative proportion of each protonation state in comparison to other states present in solution based on the pKa constants of the ionisable groups using the Henderson-Hasselbalch equation and incorporated effects of water temperature as a pKa influencing factor. The results were visualised across the pH range ( Fig. 1) and compared between today's average sea surface pH (pH 8.1), future oceanic conditions (pH 7.7) 1 and temporary coastal scenarios (pH 7.2) 5 ( Table 1). The protonated toxic form of tetrodotoxin will increase by 6.2% under the RCP8.5 scenario (pH 7.7, sea surface temperature (SST) 20.1°C), whilst the presence of the most toxic saxitoxin state, with protonated 1,2,3-and 7,8,9-guanidinium groups, will increase by 15.5% (Table 1). Taking salt (KCl) into account alters the change towards the fully protonated form of STX to 13.0% (for details see Methods). Protonated and therefore toxic forms of saxitoxin derivatives, for example neoSTX and dcSTX, also increase by 9% and 17% (see Table 1).
To investigate the electrostatic properties of the protonated toxin forms causing enhanced toxicity, we computed lowest energy models of current and future TTX and STX protonation states using our recently developed and experimentally verified quantum chemical approach 26 . We calculated the molecular electrostatic potential (MEP) without and with the presence of water molecules (see also Supplementary Information). We visualised the charge distribution of each conformer using their MEP mapped on an electron density iso-surface to highlight molecular differences. The three-dimensional conformation of TTX and STX does not change significantly upon protonation (root mean square deviation (RMSD) of carbon atoms between protonated and non-protonated forms is ± 0.013 Å). However, the TTX 0 and STX + protonation states show a distinct charge separation while the fully protonated states TTX + and STX 2+ are overall more positively charged (Fig. 1). The most significant changes in charge from negative to positive can be observed directly at the groups subject to protonation: the oxygen bound to C-10 in TTX and the 7,8,9-guanidinium group in STX, as well as at the TTX 7,8,9-guanidinium group. Addition of explicit water molecules around the toxins is shown to have no significant impact on the charge distribution pattern (see Supplementary Fig. S1). The increased positive charge at the imidazole guanidinium groups observed in our fully protonated models of both STX and TTX, matches with the proposed mechanism of enhanced molecular toxicity 24 . The increased relative proportion of active toxin, combined with a slower degradation rate of TTX in lower pH conditions 27 and minimal pH-effects on the receptors in the pH range of ocean acidification 28 , suggests a significantly increased bioavailability of these keystone molecules, and thus a significant increase of their toxicity, in future oceans.
To visualise the increased abundance of protonated toxic forms of STX in the ocean at a global scale we produced geospatial maps for current (Fig. 2a) and future oceanic pH and sea surface temperature conditions The future increase of active toxin forms predicted by our geospatial interpolation models ( Fig. 2) is relative to the total amount of toxin present. Combining this proportional increase in toxicity with the projected increases in HAB duration, intensity 17 and actual higher toxin production within the cells 33 could result in devastating effects on marine fisheries, tourism, coastal ecosystems, and public health 21 , especially once the toxins are released during HAB termination. The increase in protonated and therefore toxic forms also extends to a multitude of STX derivatives produced by HAB forming algae that vary with local environmental conditions 34 , as all main STX derivatives share the 7,8,9-guanidinium group and therefore increase in their toxicity like STX does (see Table 1 and Supplementary Table 1).
Recent years have seen rising numbers of STX-related PSP recordings from cold northern waters 19,35 , such as the Barents Sea where the STX producing Alexandrium tamarense occurs 36 . In these areas, algal toxins, in particular STX, were identified in ten out of 13 marine stranded or harvested mammal species 35 , including humpback and bowhead whales, seals and sea otters. Since many of these affected mammals prey upon filter feeders such as the Alaskan butter clam (Saxidomus gigantea), a species also frequently consumed by the local people, an increase in toxicity as indicated by our maps for this region would have even more devastating direct implications.
We therefore applied our model to calculate the projected toxicity at the end of this century using current STX contents determined in butter clams collected from an affected area and found that the amount of toxic STX in butter clams from Alaska will increase in the future to levels exceeding the current US Food and Drug Administration (FDA) limit, putting marine predators and food security at risk (Fig. 3). To maintain the current recommendations for seafood safety in the future (RCP8.5 conditions), the FDA limit of 80 g/100g total saxitoxin in tissue, which equals 50.4 g/100g of toxic STX form (Fig. 3), will need to be reduced by over 20% to 62.8 g/100g of total saxitoxin. Despite seasonal variability with clear saxitoxin summer peaks (see Supplementary Fig. S3) all butter clam samples taken since May 2014 already exceed the current FDA limit (Fig. 3b). In combination with projected increases of total STX concentration released by HABs in future ocean conditions 32 , our estimates made here for future STX toxicity may even be exceeded.
The most-impacted areas also encompass the Great Barrier Reef and the Solomon Islands (Fig 2c), where organisms ranging from dinoflagellates to worms, blue ring octopuses and pufferfish use STX and TTX as signalling molecules for key ecological functions 14 . Both toxins play a vital role in species interactions, such as deterring potential predators, attracting potential mates or serving as venom to overcome larger prey 14 . An imbalance of these interactions caused by altered effectiveness of these signalling molecules could significantly impact the ecological network. Many other biomolecules used by marine organisms to communicate have pH-sensitive chemical groups similar to the guanidinium groups shown here for TTX and STX 9 and are likely to be altered by future climate change. This impact of pH can further be expected to apply not only to biomolecules but virtually any molecule dissolved in the sea that can be protonated, from marine drugs to man-made pollutants, such as pharmaceuticals, pesticides or plasticisers. However, the responses of marine organisms under future ocean conditions can be variable 8 , and difficult to predict owing to species specific differences and their potential for adaptation. A better understanding of the impacts of pH and temperature on chemicals used by marine organisms is urgently needed to assess the full risk for marine life in changing oceans.

Calculation of protonation state abundance
Different protonation states of a molecule are present at different pH conditions. The pH at which 50% of a given ionisable group are protonated and 50% remain unchanged is given by a group-specific pKa value, which can be determined by potentiometric or NMR-based titration. 37 Figure 1. Temperature was incorporated as a factor influencing the pKa constants by -0.02 units for +1°C. This factor has been established for primary amine groups, similar to the guanidinium groups in TTX and STX, by Reijenga et al. 42 . It was used here to calculate the data plotted for the ±10°C curves framing the abundance curves calculated at the respective titration temperatures (25°C for TTX, 20°C for STX) in Fig. 1, because experimental pKa values determined at other temperatures are not available. The changes in abundance of the protonated, toxic forms in Table   1 were calculated based on the same pKa values. Differences between current and future/coastal pH where temperature changes were taken into account by employing the above described pKa-influencing temperature factor. For the RCP8.5 scenario we assumed global annual average ocean sea surface temperature (16.1°C) 44 and average pH 8.1 for current conditions and pH 7.7 (average)/ pH 7.2 (coastal) combined with a 4°C SST increase to average 20.1°C for future conditions (RCP8.5). In coastal areas pH 7.2 is already observed temporarily 45 and is likely to become more frequent in the future.

Optimisation of protonation state conformers
A change in protonation states of these molecules could be accompanied by structural changes to the cues in the lowered pH of future oceans. To investigate this, we used quantum chemical calculations to obtain the energetically most favourable conformers for each possible protonation state. These model conformers were then used to assess conformational differences between the protonation states, as well as differences in their molecular electrostatic potential (MEP), which describes the charge distribution around the molecule. Starting from the plain structure SMILES code of TTX and STX (PubChem 46 CID=11174599 47 and CID=37165 48 ), protons were added/removed according to the protonation state structures proposed by Mosher 49 for TTX (with the dissociating proton located at the hydroxyl group of C-10, not at the 7,8,9guanidinium group) and Shimizu et al. 40 for STX. Then conformers were optimized using the PBE0 exchange correlation functional 50 with a pc-2 basis set [51][52][53] and water as implicit solvent using COSMO 54 implemented in the ORCA suite of programs 55 (Version 3.0.0). We used the RIJ-COSX approximation 56 with a def2-TZVPP/J auxiliary basis set 57 and included D3 dispersion corrections following Grimme et al. 58 The VeryTightSCF and TightOpt criteria implemented in ORCA were used to stop the SCF gradient and the optimization at a total energy change of <10 -8 Eh, respectively. Differences of conformers between protonation states were assessed by calculation of the root-mean square deviation (RMSD) of atom coordinates after normalisation with respect to the position of C1.

Calculation and visualisation of the charge distribution
The calculation of the molecular electrostatic potential (MEP) was performed with the GAMESS program (vJan122009R1) using the PBE0 exchange correlation functional 50 in conjunction with a pc-2 basis set [51][52][53] . A three-dimensional electron density iso-surface was visualized with 100 grid points, a medium grid size and a contour value of 0.03 e×a0 -3 using the wxMacMolPlt program 59 (v7.5141). The density iso-surface was coloured according to the MEP with a maximum value of 0.25 Ehe -1 and the RGB colour scheme with red representing positive, green neutral and blue negative charge.

Interpolation maps for spatial prediction of protonation state
In order to visualize the spatial distribution of the current protonation levels of saxitoxin as well as the effect of future changing oceanic pH and predicted increase of sea surface temperature on these, we generated Kernel interpolation maps with standard error for current and future predicted protonation states in ArcMap (V10.5.1) based on 6485 global occurrence records of saxitoxin-related PSP HABs and HAB causing dinoflagellates Gymnodinium and Alexandrium.
These occurrence records for paralytic shellfish poisoning were obtained from the Harmful Algal Information System metadatabase (HAEDAT, http://haedat.iode.org). From this metadatabase, we selected records for HAB (Harmful Algal Blooms) involving PSP (Paralytic Shellfish Poisoning) and filtered these records for proven presence of saxitoxin. We obtained 138 unique georeferenced records for localities with STX-related HABs. Additionally, we obtained 6347 records for the distribution of two dinoflagellate genera, which are known to produce STX, Gymnodinium and Alexandrium, from the NOAA COPEPOD database. The data has been generated from 1954-1999 scientific plankton sampling expeditions. These records were likewise curated and verified by hand. Together, this gives a first, coarse estimate for the global distribution of STX-producing marine dinoflagellates.
In order to visualise future changes in estimated toxicity, we obtained raster data of current and future pH (measured as the acidity of the ocean surface), and current and future mean sea surface temperature SSH (measured as the water temperature at the ocean surface within the topmost meter of the water column in °C). Current pH and current as well as future (2100) SST were obtained from GMED (V1.0 60 ). Within this data set, the current pH layer was sourced from 1961-2009 WOD in situ measurements based on Surface Ocean Station Data (OSD), and High-resolution Conductivity-Temperature-Depth (CTD) 61 64 .
To model spatial relationships between layers, the Geostatistical Analyst toolkit was used in ArcMap. As the data is modelled within the ocean, a world vector shorelines shapefile (GSHHS_c_L1 containing all continents except Antarctica, crude shoreline) was obtained from NOAA to serve as a barrier feature. Though there were no HAB or dinoflagellate records at these six coordinates, they were located in an area of high current pH turnover and thus proved helpful to yield a more reliable interpolation model. Their location is also indicated with high model uncertainty, reflecting their substituted nature.

Calculation of future toxicity of STX in clam tissue and FDA limit
The saxitoxin content in g/100 g clam tissue was extracted from the PSP Program website of the Quagan Tayagungin Tribe 65 for the time frame between June 2012 and July 2018 for each month and averaged annually. The proportions of toxic STX form at current conditions (pH, T) as well as the future RCP8.5 scenario were extracted from the interpolation maps for the closest location to Spit Beach, Sand Point (Alaska), respectively. It was assumed that internal clam pH was close to the environmental pH due to the limited ability of bivalves to regulate their internal pH [66][67][68] . The proportional data was then used to calculate the amount of toxic STX in g/100 g clam tissue today and assuming the two future scenarios. It illustrates how the content of toxic STX in shellfish would be affected by future conditions (Fig. 3). We further calculated the amount of toxic STX currently present at the limit of 80 g/100 g clam tissue set by the US Food and Drug Administration (FDA) 69 , which is seen as safe to consume, and included it in Fig. 3.

Data availability
Source data for curves in Fig. 1

Code availability
The code used to calculate the proportions of different protonation states is available on request from the corresponding author.