The relationship here is one that comes out in geometric calculus. There, we consider the 2d plane as having 3 kinds of objects: scalars, vectors, and bivectors, which serve the role of the imaginary of complex analysis.
Here's how this works: we have a product of vectors such that $e_x e_y = -e_y e_x$, like the cross product, but $e_x e_x = e_y e_y = 1$, like the dot product. This "geometric product" produces a bivector, which represents an oriented plane, just as vectors represent oriented lines or directions.
Call $i = e_x e_y$, and let $f = u + iv$, in analogy to complex analysis. The condition $\nabla f= 0$ reproduces the Cauchy-Riemann condition. Explicitly, this is
$$\begin{align*} \nabla f &= (e_x \partial_x + e_y \partial_y )f \\ &= e_x \partial_x u + e_x e_x e_y \partial_x v + e_y \partial_y u + e_y e_x e_y \partial_y v \\ &= e_x \partial_x u + e_y \partial_x v + e_y \partial_y u - e_x \partial_y v \\ &= e_x (\partial_x u - \partial_y v) + e_y (\partial_y u + \partial_x v)\end{align*}$$
This is a well-known result: that $\nabla$ corresponds more to $\partial/\partial \bar z$ than it does with $\partial/\partial z$. What's interesting is that the integrability condition of $\nabla f = 0$ applies in higher dimensions also, even when the notion of complex differentiation is no longer clear.
To find the relationship between $f$ and a vector field $F$, multiply by $e_x$ on the right so $F = f e_x = u e_x - v e_y = U e_x + V e_y$, with $V = -v$. This generates the integrability condition
$$\nabla F = (\partial_x U + \partial_y V) + i (\partial_x V - \partial_y U)$$
The scalar part is the divergence and the bivector part is the curl. Both of these must be zero to satisfy the same integrability condition as the Cauchy-Riemann condition.
Edit: you might find that the vector field has a negative $y$-component compared to the complex function, but this is also a known result when converting between, say, velocity fields in 2d and complex functions. Geometric calculus only helps prove it mathematically.
This would be basically an application of the Lagrange inversion theorem, which states the following:
Suppose $w$ is defined as a function of $z$ by an equation of the form
$$f(z)=w,$$
where $f$ is analytic at a point $z_0$ and $f'(z_0)\neq 0$. Then it is possible to invert (or solve) the equation for $z$:
$$z=g(w),$$
on a neighborhood of $f(z_0)=w_0$, where $g$ is analytic at the point $w_0$.
Since you have (or have proved) all the hypothesis, you can just apply the theorem.
Best Answer
The angle between two curves through a point is measured by means of the tangent vectors to the curves at the point. Then, a map preserves angles (i.e., it is conformal) if the differential of the map preserves angles. The differential of the map is the linear map defined by the Jacobian matrix. As you have said, the Jacobian corresponds to a conformal linear map, thus proving that holomorphic functions are conformal.
In the case of manifolds, a map f between two manifolds M and M’ induces a linear map between the tangent space of M at p and the tangent space of M’ at f(p), which is called the differential of f at p. The map is conformal if the differential is conformal at every point. In the present case both manifolds M and M’ are the set of complex numbers.