We present a new computation of the asymptotic gravitational wave energy fluxes emitted by a spinning particle in circular equatorial orbits about a Kerr black hole. The particle dynamics is computed in the pole-dipole approximation, solving the Mathisson-Papapetrou equations with the Tulczyjew spin-supplementary-condition. The fluxes are computed, for the first time, by solving the 2+1 Teukolsky equation in the time-domain using hyperboloidal and horizon-penetrating coordinates. Denoting by M the black hole mass and by µ the particle mass, we cover dimensionless background spins a/M = (0, ±0.9) and dimensionless particle spins −0.9 ≤ S/µ 2 ≤ +0.9. Our results span orbits of Boyer-Lindquist coordinate radii 4 ≤ r/M ≤ 30; notably, we investigate the strong-field regime, in some cases even beyond the last-stable-orbit. We compare our numerical results for the gravitational wave fluxes with the 2.5th order accurate post-Newtonian (PN) prediction obtained analytically by Tanaka et al. [Phys. Rev. D 54, 3762]: we find an unambiguos trend of the PN-prediction towards the numerical results when r is large. At r/M = 30 the fractional agreement between the full numerical flux, approximated as the sum over the modes m = 1, 2, 3, and the PN prediction is 0.5% in all cases tested. This is close to our fractional numerical accuracy (∼ 0.2%). For smaller radii, the agreement between the 2.5PN prediction and the numerical result progressively deteriorates, as expected. Our numerical data will be essential to develop suitably resummed expressions of PN-analytical fluxes in order to improve their accuracy in the strong-field regime.