The efficiency of earthquake clustering investigation is improved as we gain access to larger datasets due to the increase of earthquake detectability. We aim to demonstrate the robustness of a new clustering method, MAP-DBSCAN, and to present a comprehensive analysis of the clustering properties in three major seismic zones of Greece during 2012–2019. A time-dependent stochastic point model, the Markovian Arrival Process (MAP), is implemented for the detection of change-points in the seismicity rate and subsequently, a density-based clustering algorithm, DBSCAN, is used for grouping the events into spatiotemporal clusters. The two-step clustering procedure, MAP-DBSCAN, is compared with other existing methods (Gardner-Knopoff, Reasenberg, Nearest-Neighbor) on a simulated earthquake catalog and is proven highly competitive as in most cases outperforms the tested algorithms. Next, the earthquake clusters in the three areas are detected and the regional variability of their productivity rates is investigated based on the generic estimates of the Epidemic Type Aftershock Sequence (ETAS) model. The seismicity in the seismic zone of Corinth Gulf is characterized by low aftershock productivity and high background rates, indicating the dominance of swarm activity, whereas in Central Ionian Islands seismic zone where main shock-aftershock sequences dominate, the aftershock productivity rates are higher. The productivity in the seismic zone of North Aegean Sea vary significantly among clusters probably due to the co-existence of swarm activity and aftershock sequences. We believe that incorporating regional variations of the productivity into forecasting models, such as the ETAS model, it might improve operational earthquake forecasting.