Is every scalar differential operator on $(M,g)$ that commutes with isometries a polynomial of the Laplacian?

On $(\mathbb{R}^n, g_{\text{std}})$ with $\Delta$ the Laplacian, the following holds:

Fact: Every scalar differential operator $D$ that satisfies $D \circ F^* = F^* \circ D$ for all isometries $F \in \text{Isom}(\mathbb{R}^n)$ can be written as $D = P(\Delta)$ for some polynomial $P$.

My question is whether a similar statement is true on any Riemannian manifold $(M,g)$, where $\Delta$ becomes the Laplacian with respect to the metric.

As usual, this question is a refinement of a previous question of mine.

On a (pseudo)-Riemannian manifold $(M,g)$, in addition to the metric $g_{a b}$, there is a family of natural tensors, namely the Riemannian curvature $R_{abcd}$ and its covariant derivatives $R^{(k)}{}_{a b c d p_1 \dots p_k} := \nabla_{p_1}\dots\nabla_{p_k} R_{abcd}$ for any finite $k$. All these tensors can be used as the symbols of natural Riemannian differential operators. I use the abstract index notation, see more on it here and here.

"Natural" is, roughly speaking, the synonym for "commuting with diffeomeorhisms" (see $[1]$), or "... with isometries", as in $[2]$.

The Laplacian $\Delta t := g^{ab} \nabla_a \nabla_b t$ is a differential operator with the symbol $g^{a b}$, that is the inverse metric $g^{-1}$. Here $t$ may be a function or an arbitrary tensor (with the indices suppressed). This is not the only natural differential operator on (pseudo)-Riemannian manifolds with this symbol. Such operators form a family known as generalized Laplacians.

It is well known, however, that all the (scalar) natural differential operators on a (pseudo)-Riemannian manifold can be obtained as linear combinations of (complete) contractions of monomials of the form $$ g^{-1} \otimes \dots \otimes g^{-1} \{\otimes \epsilon \otimes \} \otimes R^{(l_1)} \otimes \dots R^{(l_r)} \otimes \nabla_{a_1} \dots \nabla_{a_s} \tag{*} $$ where $\{\otimes \epsilon \otimes \}$ indicates an optional factor, the Riemannian volume density (= volume form on an orientable manifold). A proof can be found in $[3]$. See also this answer for a related question.

So, for example, the operator $Ric^{a b} \nabla_a \nabla_b$ is a natural (pseudo-)Riemannian operator, that is it commutes with all the isometries of $(M,g)$.

The mentioned fact clearly implies the theorem, cited by Jesse, because the only possibility in $(*)$ to contract the metric $g^{a_ i a_j}$ with two instances $\nabla_{a_i}$ and $\nabla_{a_j}$ of the Levi-Civita connection yields a Laplacian (in the Eulcidean case commuting the derivatives does not produce curvatures).

A modern introduction (and the canonical reference) to this circle of questions is $[4]$.

REFERENCES.

M. Atiyah, R. Bott, V. K. Patodi, On the heat equation and the index theorem, Inventiones Mathematicae, 1973, 19(4), pp. 279-330
D.B.A. Epstein, Natural tensors on Riemannian manifolds, J. Differential Geom., 1975, 10(40), pp. 631-645
P. Stredder, Natural differential operators on Riemannian manifolds and representations of the orthogonal and special orthogonal groups, J. Differential Geom., 1975, 10(4), pp. 647-660
I. Kolar, P.W. Michor, J. Slovak, Natural Operations in Differential Geometry, Springer, 1993.