Let me comment on what I know in an open set $\Omega$, trying to control $C$. First of all, the heat semigroup can be expressed through a kernel $p$ which is pointwise dominated by the heat kernel in $\mathbb{R}^n$, written by Bazin, but this is not sufficient, because of $\nabla_x$. One can extend the Gaussian bound to complex times thus getting
$$|p(z,x,y)| \le C(\Re z)^{-d/2} \exp (-c\frac{|x-y|^2}{|z|})$$ if $z=x+iy$ and, say, $|y| \le x$, with constants $C,c$ depending only on the constants in the real estimate, hence independent of $\Omega$. See Chapter 6 of the book "Analysis of Heat Equations On Domains", by El Maati Ouhabaz. The constants can be made explicit by following the proof. Next, using Cauchy theorem for holomorphic functions one deduces the estimate $$|p_t(t,x,y)| \le \frac{C}{t}p(x,y,ct),$$ see again the same book. Since $p_t=\Delta_x p$, this gives $$\|\Delta e^{t\Delta}\|_p \le \frac{C}{t}$$ with $C$ still independent of $\Omega$ and computable. The last point needs the companion "elliptic" interpolation inequality for the gradient. In the whole space it is
$$\|\nabla u\|^2_p \le C \|u\|_p\|\Delta u\|_p$$ and is usually proved using the fundamental solution of the heat equation! When $p=2$ the above inequality holds in $\Omega$, with $C=1$, by integration by parts and, with some effort, can be extended to $1\le p<2$ without requiring any assumption on $\Omega$. However, if $p>2$ I do not see how to avoid some regularity of the domain. Note that the elliptic gradient estimate is equivalent to $$\|\nabla u\|_p \le \epsilon \|\Delta u\|_p + \frac{C}{\epsilon} \|u\|_p$$ for every $\epsilon >0$. This inequality holds in $\Omega$, if $\Omega$ has the extension property and the Laplacian is substituted by the full $W^{2,p}$-norm. Then one needs elliptic regularity up to the boundary and obtains the gradient estimate in the above form for small $\epsilon$. Large $\epsilon$ are then deduced by the semigroup law, since the heat semigroup decays exponentially in $\Omega$. It is also possible to prove pointwise bounds on $\nabla_x p$ (with $1/\sqrt t$ in front) by interpolating the estimates for $p$ and $\Delta p$ in small balls, but again one needs some regularity of the boundary to do it and the control on the constants is not very explicit.