On equation $\Delta \circ \partial/\partial X=\partial/\partial X \circ \Delta$ on a Riemannian manifold
This is an expansion on Igor's comment and fixes my previous mistake (see also Willie's answer). A direct computation gives $$ D(f) = 2f^{ab}\nabla_{(a}X_{b)} + X_a(\nabla^b\nabla_b\nabla^a-\nabla^a\nabla_b\nabla^b)f + f^a\nabla^b\nabla_bX_a . $$ The middle summand is the same as $R_{ab}f^a$, while the formula $\nabla^b\nabla_aX_b=\nabla_a\nabla^bX_b+R_a{}^bX_b$ yields $$ D(f) = 2f^{ab}\nabla_{(a}X_{b)} + (\nabla^b\nabla_bX_a + \nabla^b\nabla_aX_b - \nabla_a\nabla^bX_b)f^a . $$ This can be rewritten $$ D(f) = \langle L_Xg, \nabla^2f \rangle + \langle \nabla f, \delta L_Xg - d\delta X\rangle , $$ which easily gives the formula for $D$ when $X$ is a conformal Killing field. The interpretations for your questions are as follows:
Case 1: $\nabla_{(a}X_{b)}$ is not identically zero. This is when $D$ is a second-order operator. Note that the trace of $\nabla_{(a}X_{b)}$ is the divergence. By the divergence theorem, the integral of the trace is zero, so the bilinear form $\nabla_{(a}X_{b)}$ cannot be positive definite. Thus there is no example of the type requested in your second question. (Note that on noncompact manifolds there are examples: e.g. Euclidean space with the Euler vector field $X=\sum x^i\partial_{x^i}$, so $D$ is proportional to the Laplacian.)
Case 2: $\nabla_{(a}X_{b)}\equiv0$. (That is, $X$ is Killing.) Since the trace of $L_Xg$ is $2\delta X$, we see that $D\equiv0$. That is, the Lie algebra of vector fields satisfying ($\ast$) is the Lie algebra of Killing vector fields.
The following formula is known among the experts but hard to find in the literature, so I figure I will document it here. Throughout $(M,g)$ denote an arbitrary pseudo-Riemannian manifold, and $\nabla$ its Levi-Civita connection.
Definition Given a vector field $X$, its corresponding 0th order deformation tensor is defined to be ${}^{(X,0)}\pi := \mathcal{L}_X g$, where $\mathcal{L}_X$ is Lie differentiation with respect to $X$.
The corresponding 1st order deformation tensor is defined using a formula similar to that of Christoffel symbols: $$ {}^{(X,1)}\pi_{ab}{}^c := \frac12 g^{cd} \left[ \nabla_a ( {}^{(X,0)}\pi_{bd}) + \nabla_b ({}^{(X,0)}\pi_{ad}) - \nabla_d ({}^{(X,0)}\pi_{ab}) \right] $$Lemma Let $\Xi$ be an arbitrary $k$-covariant tensor field. And let $X$ be a vector field. The following formula holds for the commutation: $$ [ \nabla_a, \mathcal{L}_X ] \Xi_{b_1\cdots b_k} = \sum_{j = 1}^k {}^{(X,1)}\pi_{a b_j}{}^c \Xi_{b_1 \cdots b_{j-1} c b_{j+1} \cdots b_k} $$
With the aid of these formulas, we have immediately that, writing $\triangle_g$ for the Laplace-Beltrami operator, first
$$ [ \nabla_X, \triangle_g] f = [\mathcal{L}_X, \triangle_g ] f $$
because Lie derivation and covariant differentiation act identically on scalars, and then
$$ [\mathcal{L}_X, g^{ab}\nabla_a\nabla_b] f = \mathcal{L}_X (g^{ab}) \nabla^a \nabla_b f + g^{ab} [\mathcal{L}_X, \nabla_a] \nabla_b f + g^{ab} \nabla_a [\mathcal{L}_X, \nabla_b ]f $$
The first factor we can compute to get
$$ \mathcal{L}_X(g^{ab}) = - {}^{(X,0)}\pi^{ab} $$
using that $g^{ab} g_{bc} = \delta^a_c$. The third factor vanishes because Lie differentiation commutes with exterior differentiation. And we use our Lemma for the second term. We get, finally
$$ [\nabla_X, \triangle_g] f = - {}^{(X,0)}\pi^{ab} \nabla_a\nabla_b f - g^{ab} ~{}^{(X,1)}\pi_{ab}{}^c \nabla_c f. $$
Remarks:
Notice that the first order deformation tensor is defined in terms of the 0th order one. So that when $X$ is Killing, automatically both the ${}^{(X,0)}\pi$ and ${}^{(X,1)}\pi$ vanish, and differentiation with $X$ commutes with the Laplacian.
In the case ${}^{(X,0)}\pi = \phi g$ for some scalar function $\phi$ (so $X$ is conformally Killing), one can check that the formula reduces to the one I gave in a comment above.
When the function $\phi$ in the previous item is a non-zero constant (which some people refer to as $X$ being a homothetic vector field) one gets the special case $$ [ \nabla_X, \triangle_g] f = \phi \triangle_g f $$