My understanding is that, in complex domain, if $f$ is a holomorphic function around $z_0$, it can be locally expressed as its own Taylor series, therefore $f(z_0+\epsilon i) \approx f(z_0)+\epsilon i f'(z_0)$. By setting $z_0=x$ and $z=i$, I can conclude the approximation asked in this question.
However, I am in a bit of doubt, especially from the geometrical point of view. What does it mean that the rate of change along the complex direction ($i$) is determined by the derivative of $f$ along the real line direction? (probably I am missing something here)
Additionally, if the above approximation correct, does $f$ has to be holomorphic for $f(x+\epsilon i) \approx f(x)+\epsilon i \frac{d}{dx}f(x)$ to be correct?
Here is the definition of the derivative of a function on the complex plane:
$$f'(z)=\lim_{|\epsilon|\to 0} \frac{f(z+\epsilon)-f(z)}{\epsilon}.$$
Here, $\epsilon$ can approach zero along any direction. This means $\epsilon$ could go along small real numbers or small imaginary numbers. A function is complex differentiable at $z$ if this limit exists. The hypothesis of complex differentiability is that the real limit and imaginary limit both exist and are well defined. This compatibility gives the famous Cauchy-Riemann equation
$$\frac{\partial u}{\partial x}=\frac{\partial v}{\partial y},\,\, \frac{\partial u}{\partial y}=-\frac{\partial v}{\partial x},$$
where $f(x+iy)=u(x,y)+v(x,y)i$. A function on the real numbers being differentiable isn't that strong of a condition, but being complex differentiable is huge. In particular, being complex differentiable in a neighborhood of $z_0$ means that $f(z)$ can be expressed as a power series in $z-z_0$.
Getting back to your approximation, the definition of $f'(z)$ gives
$$f(z+i\epsilon)=f(z)+\epsilon if'(z).$$
Since $f'(z)$ can be computed as the derivative with respect to $x$, this approximation is perfectly valid. Additionally, being holomorphic (which is defined to mean complex differentiable) is required. If it were only real differentiable, then the Cauchy-Riemann equations wouldn't hold and your concerns about the real line determining imaginary information would be valid.