Is the space of harmonic functions invariant under the derivational operator associated with a geodesible flow?

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

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

The question is related to the following post

Vector Fields in a Riemannian Manifold

• You mean $Vf = 0$? – Arctic Char Aug 11 '17 at 12:49
• No, that is not true, but likely you made a mistake in formulation --- since you mention $\Omega$ , you probably wanted to use it... – Anton Petrunin Aug 11 '17 at 18:15
• @AntonPetrunin $\Omega$ is used to define the Laplacian $\Delta=Div \circ \nabla$, so I think it is necessary and depend on the metric. but what is a counter example? – Ali Taghavi Aug 11 '17 at 18:17
• @JohnMa No I do not mean "$Vf=0"$ – Ali Taghavi Aug 11 '17 at 18:18

Take a warped product $M=\mathbb S^1\times_f\mathbb S^1$ for a nonconstant smooth function $f$.
The horizontal vector field $V$ satisfies your condition, but for harmonic function is not invariant $f$, the function $Vf$ is not harmonic.
To see this, consider the $V$-flow $\Phi^t\colon M\to M$. If $Vf$ is always harmonic then so is $f\circ\Phi^t$ for any $t$. The latter fails evidently.