Vector field with holomorphic flow
As Ben's argument suggests, the proof that, if the flow of $X$ preserves $J$ then the flow of $JX$ preserves $J$ does depend on the integrability of $J$.
As a concrete example of an almost-complex manifold for which this property does not hold, consider $S^6 = \mathrm{G}_2/\mathrm{SU}(3)$ with its standard (non-integrable) $\mathrm{G}_2$-invariant almost-complex structure $J$.
It is known (Submanifolds and special structures on the octonians, J. Differential Geom. 17 (1982), 185–232) that the vector fields on $S^6$ whose flows preserve $J$ consist exactly of the vector fields associated to $1$-parameter subgroups of $\mathrm{G}_2$, and hence they form a $14$-dimensional Lie algebra ${\frak{g}}_2$. Moreover, any vector field $X$ on a connected open subset $U\subset S^6$ whose flow preserves $J$ is the restriction to $U$ of a (unique) global vector field on $S^6$ that preserves $J$.
However, any such global $X$ that is not identically zero must vanish at some point $p\in S^6$, and its linearization at $p$ is a nonzero linear transformation $X'(p):T_pS^6\to T_pS^6$ with purely imaginary eigenvalues (since the flow of $X$ must also preserve the metric on $S^6$). It follows that the linearization of $JX$ at $p$ is a nonzero linear transformation with real eigenvalues. Consequently, $JX$ cannot belong to ${\frak{g}}_2$.
A vector field $X$ has flow preserving $J$ just when $L_X J=0$, i.e. just when $[X,JY]=J[X,Y]$ for any vector field $Y$. Take $Y$ to also have such flow. Many such $Y$ exist locally, spanning the tangent bundle, because $(M,J)$ is a complex manifold. Vanishing of the Nijenhuis tensor gives $[JX,JY]=[X,Y]+J([JX,Y]+[X,JY])$. But then we apply the flows of $X$ and $Y$: $[JX,JY]=[X,Y]+J(J[X,Y]+J[X,Y])=-[X,Y]=J^2[X,Y]=J[X,JY]=J[JX,Y]$. So $[JX,JY]=J[JX,Y]$. So $L_{JX} J=0$.
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?