Why is a diffeomorphism an isometry if and only if it commutes with the Laplacian?

I came across the following statement in a book on automorphic forms:

In general, on a Riemannian manifold, the Laplace-Beltrami operator $\Delta$ is characterized by the property that a diffeomorphism is an isometry if and only if it leaves $\Delta$ invariant.

For this question, I am interested in the last part of the statement, that "a diffeomorphism is an isometry if and only if it leaves $\Delta$ invariant." What is a precise statement of this result? We're viewing $\Delta$ as an operator on $C^\infty(M)$, right? Is there a way to prove it just from the description of $\Delta$ in coordinates: $$\frac{-1}{\sqrt{\det{g}}}\sum_{i,j}\partial_j(g^{ij}\sqrt{\det{g}}\partial_i f)?$$ I'd also be interested in a way to prove it from the definition of $\Delta$ as $-\text{div}\circ\nabla$.

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.