To do De Giorgi-Nash-Moser it is important that $f$ is in $L^q$ for $q$ larger than half the dimension. A heuristic is scaling: if $v$ solves $-\Delta v = f$, then the right side for the rescaling $v(\epsilon x)$ is $\epsilon^2 f(\epsilon x)$, whose $L^q$ norm is like $\epsilon^{2-n/q}$, so zooming in helps when $q > n/2$.

When in dimension $4$ or higher, one can take e.g. $v = \sum h_k\varphi(2^k x)$ for some smooth $\varphi$ that is $1$ in $B_1$ and vanishes outside $B_2$. Taking $h_k = 1/k$ makes $v$ unbounded ($\log\log$ growth near the origin) while $f := -\Delta v$ is $L^2$ (the integral of $f^2$ is like $\sum h_k^2 2^{(4-n)k}$).

It seems to me that for such a choice of $f$, the solution $u$ will immediately become unbounded since $u - v$ is unbounded initially and solves the heat equation.

Note also that indeed, in Deane's computation, the (reciprocal) exponent on $f$ looks like $n/2$ for large $p$. It is important to use that $f \in L^q$ for $q > n/2$, for instance as follows (in the elliptic case for simplicity): Let $A_p = \left(\int u^{p\chi}\right)^{\frac{1}{p \chi}}$, with $\chi = n/(n-2)$. Applying Holder to the last term in the first line of Deane's computation gives an inequality like
$$A_p \leq (pC(f,S))^{1/p} A_{\gamma p},$$
where $\gamma < 1$ depends on $q/(q-1) < \chi$, and $C(f,S)$ depends on $f$ and the Sobolev constant. Iteration gives
$$A_{2\gamma^{-k}} \leq C(f,S)^{\sum j\gamma^j} \|u\|_{L^2(Q)}.$$
Taking $k \rightarrow \infty$ gives an $L^{\infty}$ bound in terms of the desired quantities and the Sobolev constant. I'm not sure how to remove dependence on the Sobolev constant.