Harmful algal blooms (HABs) negatively impact coastal ecosystems, fisheries, and human health, and their prediction has become imperative for effective coastal management. This study aimed to evaluate spatial-temporal variability patterns and phenology for key toxigenic phytoplankton species off southern Portugal, during a 6-year period, and identify region-specific environmental drivers and predictors. Total abundance of species responsible for amnesic shellfish poisoning (Pseudo-nitzschia spp.), diarrhetic shellfish poisoning (Dinophysis spp.), and paralytic shellfish poisoning (G. catenatum) were retrieved, from the National Bivalve Mollusk Monitoring System public database. Contemporaneous environmental variables were acquired from satellite remote sensing, model-derived data, and in situ observations, and generalized additive models (GAMs) were used to explore the functional relationships between HABs and environmental variables and identify region-specific predictors. Pseudo-nitzschia spp. showed a bimodal annual cycle for most coastal production areas, with spring and summer maxima, reflecting the increase in light intensity during the mixed layer shoaling stage, and the later stimulatory effects of upwelling events, with a higher bloom frequency over coastal areas subjected to stronger upwelling intensity. Dinophysis spp. exhibited a unimodal annual cycle, with spring/summer maxima associated with stratified conditions, that typically promote dinoflagellates. Dinophysis spp. blooms were delayed with respect to Pseudo-nitzschia spp. spring blooms, and followed by Pseudo-nitzschia spp. summer blooms, probably reflecting upwelling-relaxation cycles. G. catenatum occurred occasionally, namely in areas more influenced by river discharges, under weaker upwelling. Statistical-empirical models (GAMs) explained 7-8%, and 21-54% of the variability in Pseudo-nitzschia spp. and Dinophysis spp., respectively. Overall, a set of four easily accessible environmental variables, surface photosynthetically available radiation, mixed layer depth, sea surface temperature, and chlorophyll-a concentration, emerged as the most influential predictors. Additionally, over the coastal production areas along the south coast, river discharges exerted minor negative effects on both HAB groups. Despite evidence supporting the role of upwelling intensity as an environmental driver of Pseudo-nitzschia spp., it was not identified as a relevant model predictor. Future model developments, such as the inclusion of additional environmental variables, and the implementation of species- and period-specific, and hybrid modelling approaches, may further support HAB operational forecasting and managing over complex coastal domains. Copyright © 2022 Elsevier B.V. All rights reserved.