Assume that $V$ is a vector field on a Riemannian manifold $(M,g)$ with natural volum form $\Omega$ arising from $g$. Assume that the solution curves of $V$ are parametrized geodesics of the Riemannian meteic $g$.

Is it true to say that the space of harmonic functions are invariant under the derivational operator $D(f)=V.f=df(V)$?