Understanding the spatial variability of forest gross primary productivity (GPP) is essential for diagnosing ecosystem functioning and supporting monitoring and management within protected areas. However, many data-driven approaches emphasize predictive performance while offering limited interpretability of the drivers underlying spatial heterogeneity. Here, we develop an interpretable machine-learning framework to diagnose spatial patterns and dominant drivers of forest GPP within China's national-level protected areas. Satellite-derived GPP and environmental variables were aggregated to a 0.1 degrees spatial grid over 1990-2018 to characterize long-term mean forest productivity. Multiple machine-learning models were evaluated, and the best-performing model was interpreted using explainable artificial intelligence to quantify driver importance and response behavior. The mean forest GPP was 759.5 g C m-2 yr-1, with pronounced spatial heterogeneity. Approximately 22% of protected areas experienced increases in forest GPP exceeding 20% between 1990 and 2018, primarily in humid regions, whereas 12.5% showed declines greater than 20%, mainly in the southern Qinghai-Tibet Plateau, northern arid regions, and west-central temperate semihumid zones. Among the evaluated models, XGBoost achieved the highest predictive performance on independent test data (R2 = 0.76, RMSE = 262 g C m-2 yr-1). Precipitation, temperature, and solar radiation emerged as the dominant drivers, with precipitation explaining 53.4% of the study area, followed by temperature (19.7%) and solar radiation (16.0%). Forest fragmentation exhibited a predominantly negative association with forest GPP. This study is a spatial diagnostic analysis that provides transparent insights to support spatial prioritization, monitoring design, and management planning within protected areas.