Let $\phi: M \to N$ and suppose $(M,g)$ and $(N,h)$ are Riemannian manifolds. If $\phi^* \circ \Delta_h = \Delta_g \circ \phi^*$, then, as Jack Lee indicated, $\phi$ is an isometry. Here is an elementary and complete proof using the $-div \circ \nabla$ point of view:
By hypothesis, for each pair of smooth functions
$u,v: N \to {\mathbb R}$ we have
\begin{eqnarray*}
%
\int_M (\Delta_g (u \circ \phi)) \cdot (v \circ \phi) \, dV_g
&=&
\int_M \left((\Delta_h u) \circ \phi \right) \cdot (v \circ \phi) \, dV_g \\
&=&
\int_N (\Delta_h u) \cdot v \cdot (\phi^{-1})^{*}(dV_g) \\
&=&
\int_N (\Delta_h u) \cdot v \cdot \rho\, dV_h \\
%
\end{eqnarray*}
where $(\phi^{-1})^{*}(dV_g) = \rho \cdot dV_h$.
If we assume that $u$ is compactly supported, then by integrating
each side by parts we obtain
$$
\int_M \bar{g} (d( u \circ \phi), d(v \circ \phi) ) \, dV_g~
=~
\int_N \bar{h} (d u, d(v \cdot \rho) ) \, dV_h.
$$
where $\bar{g}$ and $\bar{h}$ are the co-metrics associated to $g$ and $h$ respectively.
By changing variables in the integral on the left hand side, we obtain
$$
\int_N (\phi^{-1})^*(\bar{g}) (du, dv )\cdot \rho \cdot dV_h~
=~
\int_N \bar{h} (d u, d(v \cdot \rho) ) \, dV_h.
$$
If we choose $v=1$, then the left hand side equals zero and so
$$
0~
=~
\int_N \bar{h} (d u, d \rho ) \, dV_h.
$$
Since $u$ is an arbitrary smooth (compactly supported) function, this implies
that $d \rho=0$ and so $\rho$ is constant.
Therefore we may divide both sides of the previous identity by
$\rho$ and obtain
$$
\int_N (\phi^{-1})^*(\bar{g}) (du, dv )\cdot dV_h~
=~
\int_N \bar{h} (d u, dv )\cdot dV_h.
$$
It follows that $(\phi^{-1})^*(\bar{g})= \bar{h}$ and so $\phi^*(h)=g$.