I think by "standard subsolution estimate", Kurata means the following:
Let $u\ge0$ be a subsolution in $Q_2:=B_2(0)\times(-4,0)$, i.e. $\partial_t u - \nabla\cdot(A\nabla u)\le 0$. Then $u$ satisfies $$sup_{Q_{1/2}} u \le C\Big( \int_{Q_1} |u|^2 \;dxdt \Big)^{\frac{1}{2}}$$ for some constant depending only on the ellipticity contrast $\lambda$ and the dimension $d$.
Then the result of Kurata with a constant independent of $r$ follows by scaling space with $r$ and time with $r^2$. Likewise the constant may be chosen independently of $t_0$ and $x_0$ by translation invariance of the problem.
I am guessing that Kurata assumes $u\ge 0$ and $V\ge 0$, so that a solution to his problem is indeed a subsolution: $\partial_t u - \nabla\cdot(A\nabla u)\le 0$.
The above "standard subsolution estimate" may be obtained by parabolic Moser iteration, cf. J. Moser, A harnack inequality for parabolic differential equations, Comm. Pure Appl. Math., link.
Moser iteration corresponds to testing the equation with clever test functions. The simplest example would be the following: consider an elliptic problem (e.g. $u$ above does not depend on time) $-\nabla\cdot A \nabla u \le 0$ in $B_2$ and test with $\eta^2 u$, where $\eta$ is a test function which is one in $B_{1/2}$ and zero outside $B_1$. Then after integration by parts $$0\ge-\int \eta^2 u \nabla\cdot A \nabla u \;dx = \int \Big(\eta^2 \nabla u \cdot A \nabla u + 2\eta u\nabla\eta \cdot A \nabla u \Big) \;dx.$$ Ellipticity and Young's inequality yield $$\int \eta^2 |\nabla u|^2 \;dx\le C(\lambda) \int |\nabla\eta|^2 u^2 \;dx.$$ This is usually called Caccioppoli's inequality. By the product rule, it can be upgraded to $$\int |\nabla (\eta u)|^2 \;dx\le C(\lambda) \int |\nabla\eta|^2 u^2 \;dx.$$ Now we use the Sobolev embedding $\|v\|_{L^p(R^d)}\le C(d) \|\nabla v\|_{L^2(R^d)}$ with $p=2d/(d-2)$ (let's assume $d>2$) to obtain $$\Big(\int |\eta|^p |u|^p \;dx\Big)^{\frac{2}{p}}\le C(d,\lambda) \int |\nabla\eta|^2 u^2 \;dx.$$ Finally, by the choice of test function, e.g. $|\nabla \eta| \le C(d)$ and we obtain $$\Big(\int_{B_{1/2}} |u|^p \;dx\Big)^{\frac{2}{p}}\le C(d,\lambda) \int_{B_1} u^2 \;dx.$$ This was the first step to give you some gain of integrability (giving up a little bit on the domain of integration on the right hand side). For the next step, you test with something like $\eta^2u^{p-1}$...
I hid a lot of technical details and the argument gets more complicated when you iterate higher and if you do parabolic estimates, but the idea is there. Anyway, I've rambled on much longer than I intended to...
PS: The reason why the constant in Aronson and Serrin seems to depend on $V$ and $r$ is that they consider a (much) more general problem which in particular is not scaling invariant in the same way that $\partial_t u - \nabla \cdot A \nabla u \le 0$ is.