  In this paper, we prove local H\"older continuity for the spatial gradient of
weak solutions to $$u_t - \text{div} (|\nabla u|^{p-2}\nabla u) + \text{P.V.}
\int_{\mathbb{R}^n} \frac{|u(x,t) - u(y,t)|^{p-2}(u(x,t)-u(y,t))}{|x-y|^{n+ps}}
\ dy = 0.$$ It is easy to see that parabolic quasilinear equations are not
scaling invariant and this led to the development of the method of intrinsic
scaling by E.DiBenedetto, E.DiBenedetto-Y.Z.Chen, J.Kinnunen-J.Lewis and
A.Friedman-E.DiBenedetto. In a very recent paper, C.de Filippis-G.Mingione
proved gradient H\"older continuity for mixed local-nonlocal quasilinear
elliptic equations and in this paper, we extend this result to the parabolic
case. Since we only expect regularity for $\nabla_x u$ in the parabolic
setting, it is not clear how to extend the elliptic proof to the parabolic
case. In order to overcome this difficulty, we instead follow the ideas
developed by T.Kuusi-G.Mingione combined with the novel tail estimates of
C.deFilippis-G.Mingione. An advantage of our approach is that we can obtain
both $C^{1,\alpha}_x$ regularity as well as $C^{0,1} _x$ potential estimates in
one go. Moreover, we do not need to make use of any form of Caccioppoli
inequality and instead, the regularity is obtained only through a suitable
difference estimate.
