All species were projected to undergo a considerable decrease in mean habitat suitability across most of the latitudinal profile, regarding the present-day, leading to potential extirpations across (when the species is currently present) equatorial, tropical, and temperate regions. Also, despite decreasing along most of the latitudinal profile, mean habitat suitability is seen increasing over time and across RCP scenarios for all species at higher latitudes (i.e., the polar and subpolar regions). These changes in habitat suitability, then, result in specific distributional changes which, despite being different between species, share a considerable northward shift in the distributional centroid between the present-day and (at least) the middle of the century (followed by some relative oscillations for P. australis and P. fraudulenta), and the year 2100 (for P. seriata). These results hint at potential poleward distributional shifts until the end of the century for these Pseudo-nitzschia species. Concerning variable contribution, bathymetry was the most contributing variable for all species, with values between approximately 61% (P. fraudulenta) to approximately 82% (P. seriata). This trend has also been observed in SDMs with PST-producing dinoflagellates [54] and is to be expected when dealing with neritic species, and with occurrence data obtained mostly from coastal monitoring programs. At the same time, Pseudo-nitzschia blooms are known to be linked to upwelling regions [55, 56]. Afterwards, then, temperature layers were, overall, the most ecologically relevant predictors. Indeed, rising salinity has been shown to favor Pseudo-nitzschia spp. abundances [57], and as such, salinity minimum held a significant contributive role for the models of P. fraudulenta.
Rising temperatures have been linked to poleward shifts in the thermal niches of phytoplankton until the end of the century, akin to a wide variety of marine taxa [58, 59], suggesting sharp declines in tropical phytoplankton diversity in the absence of adaptation to warming [60]. Indeed, several prior studies projecting climate change biogeographical effects suggest similar phenomena could occur on phytoplankton species. Indeed, dinoflagellates in the Atlantic and Pacific oceans have been found to closely follow the rate of isotherm movement [10, 61–63], as the range of optimum environmental conditions shifts, effectively pressuring their ability to survive and adapt at lower latitudes and opening the possibility of HAB events in newer regions. Borges et al. [54] also projected potential poleward shifts for PST-producing dinoflagellate species (i.e., Alexandrium minutum; A. catenella; and Gymnodinium catenatum), while other groups such as Cocolitophores [64] also exhibit the same distributional changes. Changes in temperature have also been shown to induce significant impacts in diatoms, with warming being linked to lower cell yield and promoted growth in Pseudo-nitzschia species [65], while extreme events such as heatwaves have also triggered pelagic HAB events [37]. Indeed, a record-setting HAB event by Pseudo-nitzschia occurred as a direct consequence of a marine heatwave in the northeast Pacific Ocean between 2013–2015 [36]. Additionally, North Atlantic diatoms have been projected to shift the central positions of their core range poleward, leading to a significant reshuffling of the phytoplankton communities, with broad effects on food webs and biogeochemical cycles [30]. The present study presents similar results, particularly regarding the cold-water species P. seriata, which is projected to further restrict its southern limit of distribution in the northern hemisphere, and consequently expand its range in the Arctic and sub-Arctic regions. Despite exhibiting a similar trend, P. australis and P. fraudulenta exhibited a relatively smaller shift north, expanding more limitedly and in northern temperate and sub-arctic regions, but also exhibiting a significant restriction in lower latitudes. Diatoms are known to have a lower activation energy compared to other phytoplankton groups [66], and tropical species have been shown to adapt to multigenerational exposure to warming through various thermal strategies – i.e., either by increasing optimal growth temperature and maximum growth rate, or by shifting from specialist to generalist, increasing the maximum critical thermal limit – albeit trading off on photosynthetic efficiency, high irradiance stress, and lower growth rate [67]. Despite their potential to adapt, the lower competitive fitness of diatom species induced by thermal stress could eventually lead to species extirpation in lower latitudes, where the ocean becomes too warm to support growth [39]. The concomitant migration of Pseudo-nitzschia into new ecosystems at higher latitudes, alongside their environmental thermal optimum, poses a significant risk to marine communities and to human societies inhabiting these regions [6, 35], with these new regions supporting a widening window of optimal conditions for blooms to develop [5, 68]. However, since HAB events result from a complex interplay of abiotic and biotic interactions, the present results do not present a prediction of likelihood of ASP events, contributing instead towards the projection of species’ movements and identifying potential areas at risk.
Lastly, like any modelling endeavor, it is required to address the existence of potential caveats of the SDM analysis, which arise from methodological choices and assumptions. First, the models were subject to instances of over- and underprediction for all species. Overprediction in SDMs is common, since the models assume that species will completely occupy areas that are calculated to be climatically suitable [42], ignoring other factors – such as unaccounted environmental predictors or geographical barriers, etc. – which prevent species occurrence at a local or regional scale [69]. Concurrently, failing to predict the presence of a species in areas with existing historical records (i.e., underprediction) is directly linked to the nature of the occurrence point data. This data is retrieved mainly from online databases such as GBIF and OBIS (i.e., the Ocean and Biodiversity Information System), which do not provide a complete inventory of a species’ known distribution – mainly due to geographical biases associated with sampling efforts, due to accessibility and funding [70]. Nonetheless, instances of over or underprediction do not diminish, given that they are not dominant in the model predictions, does not diminish the predictive potential in suggesting global trends of future poleward shifts or major extirpations of species’ distribution. Second, long-term SDM projections are unable to model seasonal changes in abundance in each species, which directly impact the existence of an algal bloom. Nevertheless, the present article suggests an overall decrease in habitat suitability, followed by changes in the distribution ranges of the three species modelled into areas where they previously did not occur, or only occurred locally. These shifts pose a serious ecosystem risk, opening the possibility of nuisance blooms in novel areas [71]. Lastly, single-trophic level SDMs lack incorporating the contribution of other trophic levels to the potential expansion or restriction in distribution induced by, for instance, predator species. Indeed, phytoplankton are vulnerable to strong predation pressure by many protozoan and metazoan grazers [72], which has been shown to influence toxin production in toxic (e.g., P. seriata) and previously considered non-toxic Pseudo-nitzschia species (e.g., P. obtusa) [73, 74]. As such, including this type of trophic interactions in the models – by projecting the distribution of grazer species under the same RCP scenarios – could lead to more accurate estimates [75, 76].