Without lose of generality we may work locally in $\mathbb{R}^{2n}$ with its standard complex structure $J$. Let $\phi$ be the flow of the vector field $x'= f(x)$. So the statement in the question is a consequence of the variational equation $$ \partial_t \partial_x \phi_t = Df(\phi). \partial_x \phi_t$$
If the flow is holomorphic then $\partial_x \phi_t$ and its time derivative $\partial_t \partial _x \phi_t $ commute with $J$. This implies that $Df.J= J.Df$. Q.E.D

In the same manner one can shows the following:

**Fact:** If the flow of a vector field $x'=f(x)$ on $\mathbb{R}^n$ is harmonic then $f$ is harmonic, too.

**Proof:**

The Laplacian, the gradient and the Hessian of a vector valued function is defined naturally by componentwise corresponding operators. Hence the Hessian of a vector valued function is a $3$ dimensional matrix whose trace is a "vector".

With such notation we have the following formula:

$$\partial_t \Delta \phi= trace( {(\partial _x \phi_t)}^{tr}Hess(f)\partial_x \phi_t)+ \Delta \phi . \nabla f$$

This shows (using the time independence of $f$ to evaluate at $t = 0$) that if the flow of a vector field is harmonic then the vector field is a harmonic vector valued function.

**But what about the converse? Is the flow of a harmonic vector field, a harmonic function? Moreover, how can we rephrase the above Euclidean fact in an abstract Riemannian manifold?**