Remote estimation of phytoplankton primary production (PP) based on satellite ocean color observations has long been recognized as a key approach for investigating how marine carbon fixation responds to global climate change. The theory-based primary production model (TPM), one of the earlier proposed models, is broadly applicable across diverse aquatic ecosystems because it explicitly links PP to phytoplankton photosynthetic traits. Therefore, its accuracy is highly dependent on the performance of remote assignment of a suite of photosynthetic parameters used to calculate realized photosynthetic rates. In this study, we apply a novel and well-validated machine learning algorithm to retrieve two key photosynthetic parameters, i.e., the assimilation number (P^B_max) and the light saturation parameter (E_k) from ocean color observations. These retrievals are subsequently integrated into the TPM to generate a global, long-term monthly PP product from 1998 to 2025. Bias-correction of multi-source satellite input products and gap filling of missing data in the generated PP product are implemented to reduce uncertainties and enhance data reliability. In addition, PP products generated from alternative satellite-based models are included as benchmarks for intercomparison. Our analyses reveal contrasting spatiotemporal patterns of basin-scale PP driven by different environmental factors under climate change. The interannual variability of PP in the global ocean, as well as in the Pacific and Indian Oceans, is greatly attributed to transitions between El Niño and La Niña phases. In contrast, changes in sea-ice extent exert a dominant control on PP variability in the Arctic and the Southern Oceans. The Atlantic Ocean exhibits pronounced interannual fluctuations of PP, which may be partially explained by variability in phytoplankton community composition and chlorophyll-a concentration.