Dynamic Harmonic Regression (DHR) models are applied here to the investigation of the interannual changes in the trend and seasonality of biogeochemical variables monitored in coastal areas. A DHR model can be regarded as a time-series component model, where the phases and amplitudes of the seasonal component, as well as the trend, are parameters that vary with time, reflecting relevant changes in the evolution of the biogeochemical variables. The model parameters and their confidence bounds are estimated by data assimilation algorithms, i.e. the Kalman filter and the Fixed Interval smoother. The DHR model structure is here identified by a preliminary spectral analysis and a subsequent minimization of the Bayesian Information Criterion, thus avoiding subjective choices of the frequencies in the seasonal component. The methodology was applied to the investigation of the long-term and interannual variability of ammonia, nitrate, orthophosphate and chlorophyll-a monitored monthly in the lagoon of Venice (Italy) during the years 1986-2008. It was found that the long-term evolutions of the biogeochemical variables were characterized by non-linear patterns and by statistically significant changes in the trend. The seasonal cycles of all the variables were characterized by a marked interannual variability. In particular, the changes in the seasonality of chlorophyll and nitrate were significantly related to the changes in the seasonality of water temperature at the study site and of nutrient concentrations in river discharges, respectively. These results indicate that the methodology could be a sound alternative to more traditional approaches for investigating the impacts of changes in environmental and anthropogenic forcings on the evolution of biogeochemical variables in coastal areas.