For the maximally symmetric case with the second order operator, the proof is very simple.

A second order differential operator acting on scalars can be written in the form
$$ b^{ij} \nabla^2_{ij} + c^i \nabla_i + e $$
Consider first the action of all isometries that fixes a point $p\in M$. You must have that $b^{ij}(p)$ is a symmetric two tensor that is fixed by the action of the entire special orthogonal group, and so must be proportional to the inverse metric. Now as the principal part is the Laplacian and is invariant under the isometries, you can then conclude similarly that $c^i$ is a vector that is invariant under the action of the entire special orthogonal group, and hence vanishes.

This implies that your differential operator is
$$ \lambda \triangle_g + e $$
for scalar functions $\lambda$ and $e$. But then translation symmetry tells you immediately that both $\lambda$ and $e$ are constant.

The maximal symmetry is used crucially in the argument: for example, consider the standard torus $\mathbb{T}^2$; it is homogeneous but does not admit any rotation symmetry. In standard coordinates any second order differential operator of the form $a_1 \partial^2_{11} + b \partial^2_{12} + a_2 \partial^2_{22}$ commute with all the Killing vector fields (being spanned by $\partial_1$ and $\partial_2$).

Notice that in addition to the "natural" differential operators mentioned in the link in Phillip Andreae's comments, there are also other "less natural" differential operators that can commute with the Killing vector fields.

For example, let $(M,g)$ be any Riemannian manifold with only one Killing vector field $X$ (isometry group being $\mathbb{R}$ or $U(1)$; think a generic surface of revolution). Then quite clearly the second order operator $\mathcal{L}_X^2: f \mapsto X(X(f))$ commutes with the (one and only) Killing vector field, yet is not the Laplacian.