(Analysis is not my thing but after giving it some thought I think something like this will work — please correct me if I'm wrong)
Let $D = P(D^1, ..., D^n)$ be a differential operator commuting with isometries.
In particular $D(f \circ g) = Df \circ g$ (*) for any Schwartz function $f$ and isometry $g$.
Taking the Fourier transform $\mathfrak{F}$ we get (using the scaling property of the Fourier transform mentioned here for linear operators — I am omitting the determinant which will factor out):
$$
\mathfrak{F}(D(f \circ g)) = P \cdot \mathfrak{F}(f \circ g) = P \cdot \mathfrak{F}(f) \circ g^{-1} \\
\mathfrak{F}(Df \circ g) = F(D f) \circ g^{-1} = (P \cdot \mathfrak{F}(f)) \circ g^{-1}
$$
By (*) we get $P = P \circ g^{-1}$ (considering a $f$ such that the support of $\mathfrak{F}(f)$ is the whole real line).
The only polynomials $P(X_1, ..., X_n)$ which are invariant by all isometries are those which are functions of ${X_1}^2 + ... + {X_n}^2$ (we can show that it factors as $P - \lambda = Q \cdot ({X_1}^2 + ... + {X_n}^2)$ for some scalar $\lambda$ then the same for $Q$ etc. until we reach the degree of $P$).
Note that this should also work in the case of the Lorentz group and and the d'Alembert operator.