If $\nabla$ is the Levi-Civita connection, then they are the same (on scalar functions). This follows from the identity
$$\Gamma^i_{ij}=|g|^{-1/2}\partial_j|g|^{1/2}$$
for the contracted Christoffel symbols.
We therefore obtain
\begin{equation}
\begin{split}
g^{ij}\nabla_i\nabla_j\phi&=g^{ij}\nabla_i\partial^j\phi \\
&=\nabla_i\partial^i\phi=\partial_i\partial^i\phi+\Gamma^i_{ij}\partial^j\phi \\
&=\partial_i\partial^i\phi+\big(|g|^{-1/2}\partial_i|g|^{1/2}\big)\partial^i\phi\\
&=|g|^{-1/2}\big(\partial_i\partial^i\phi\big) |g|^{1/2} +\big(|g|^{-1/2}\partial_i|g|^{1/2}\big)\partial^i\phi\\
&=|g|^{-1/2}\partial_i\big(|g|^{1/2} \partial^i\phi\big).
\end{split}
\end{equation}
Here is a longer derivation that requires mostly basic calculus.
The Laplacian is invariant under rotations. Specifically, if $f:\mathbb R^n\rightarrow \mathbb R^n$ and $\nabla^2 f(x)\;=\alpha$, then $\nabla^2 g\;
(R^{-1} x)=\alpha$ where $R$ is any rotation matrix ($R^TR=I$ and $\det(R)=1$) and $g(x) := f(Rx)$.
Now suppose that $f:\mathbb R^n\rightarrow \mathbb R^n$ has the property that $f(x)= h(||x||)$ where $h:\mathbb [0,\infty) \rightarrow \mathbb R$, $h$ is twice differentiable, and $||x||=\sqrt{\sum_i x_i^2}.$
At $x$, we can form an orthonormal basis $\{x/||x||, v_2, v_3, \ldots v_{n-1}\}$. Now
$$
\nabla^2 f(x) = \sum_{i=1}^n \frac{\partial^2 f}{\partial x_i^2}(x)
$$
and by invariance under rotation, we can use the basis $\{x/||x||, v_2, v_3, \ldots v_{n-1}\}$ to get
\begin{equation}
(1)\quad\quad\nabla^2 f(x)= h''(||x||) + \sum_{i=2}^n p_i
\end{equation}
where $p_i:= k_i''(0)$ and $k_i(\epsilon) := f(x+\epsilon v_i)$.
The fact that $f$ only depends on the radius implies that all of the $k_i$ are the same and for $i=2, 3, \ldots n$,
$$f(x+\epsilon v_i) = k_i(\epsilon) = k(\epsilon) = f(x+ \epsilon v_2)=h(\sqrt{||x||^2+\epsilon^2}).$$
Let $r=||x||$. Using calculus,
$$k'(\epsilon) = h'(\sqrt{r^2+\epsilon^2})\frac{\epsilon}{\sqrt{r^2+\epsilon^2}}, \quad \mathrm{and} $$
$$k''(\epsilon) = h''(\sqrt{r^2+\epsilon^2})\frac{\epsilon^2}{r^2+\epsilon^2} + h'(\sqrt{r^2+\epsilon^2})\frac{\sqrt{r^2+\epsilon^2} -\epsilon \frac{\epsilon}{\sqrt{r^2+\epsilon^2}} }{r^2+\epsilon^2}.$$
Thus $p_i:= k''(0) = h''(r)\cdot 0 + h'(r)\frac{r}{r^2}=h'(r)/r$. Applying this to equation (1) gives
$$
\nabla^2f(x)= h''(||x||) + \sum_{i=2}^n p_i = h''(r) + (n-1)\frac{h'(r)}{r}.
$$
Best Answer
Here's one way to make it precise.
Theorem. Suppose $(M_1,g_1)$ and $(M_2,g_2)$ are Riemannian manifolds, $\Delta_1,\Delta_2$ are their respective Laplace operators, and $\phi\colon M_1\to M_2$ is a diffeomorphism. Then $\phi$ is an isometry if and only if for every $f\in C^\infty(M_2)$ we have $\Delta_1(f\circ\phi) = (\Delta_2 f)\circ \phi$.
(Of course, if you're only interested in diffeomorphisms from a Riemannian manifold to itself, you can take $M_1=M_2=M$ and $g_1=g_2=g$.)
I don't have time to write down a complete proof, but the idea is that the principal symbol of a differential operator is a diffeomorphism-invariant quadratic function on the cotangent bundle, and the principal symbol of the Laplace operator is given by the squared norm function on covectors. By the polarization identity, the squared norm function determines the metric on covectors (the "dual metric"), and by the musical isomorphism between $TM$ and $T^*M$, this in turn determines the metric.