Given a domain ${\rm\Omega}$ of a complete Riemannian manifold ${\mathcal{M}}$, define ${\mathcal{A}}$ to be the Laplacian with Neumann boundary condition on ${\rm\Omega}$. We prove that, under appropriate conditions, the corresponding heat kernel satisfies the Gaussian upper bound $$\begin{eqnarray}h(t,x,y)\leq \frac{C}{[V_{{\rm\Omega}}(x,\sqrt{t})V_{{\rm\Omega}}(y,\sqrt{t})]^{1/2}}\biggl(1+\frac{d^{2}(x,y)}{4t}\biggr)^{{\it\delta}}e^{-d^{2}(x,y)/4t}\quad \text{for}~t>0,~x,y\in {\rm\Omega}.\end{eqnarray}$$ Here $d$ is the geodesic distance on ${\mathcal{M}}$, $V_{{\rm\Omega}}(x,r)$ is the Riemannian volume of $B(x,r)\cap {\rm\Omega}$, where $B(x,r)$ is the geodesic ball of centre $x$ and radius $r$, and ${\it\delta}$ is a constant related to the doubling property of ${\rm\Omega}$. As a consequence we obtain analyticity of the semigroup $e^{-t{\mathcal{A}}}$ on $L^{p}({\rm\Omega})$ for all $p\in [1,\infty )$ as well as a spectral multiplier result.