Your approach is correct, I would use the mean value property to try to derive an analogue of the Cauchy estimates. Let $\Omega$ be the domain, and let $u$ be harmonic. See if you can show that if $B(z,r) \subset \Omega$, then
$$ |\partial_i u(z)| \le Cr^{-1} \|u \|_{L^{\infty}(\partial B(z,r))}$$
You will be using both the mean value property (since $\partial_i u$ is itself harmonic) and the divergence theorem. If you have further questions I can elaborate more.
Elaboration: The $L^{\infty}$, in the context of continuous functions (such as the harmonic functions in this problem), is just the norm corresponding to uniform convergence. So the sequence of harmonic functions $u_n$ converging uniformly on a compact set $K$ can be rewritten as $\|u - u_n\|_{L^{\infty}(K)} \rightarrow 0$ as $n \rightarrow \infty$.
To obtain the estimate in question, first use the mean value property to obtain:
$$ \partial_i u(z) = \frac{1}{\pi r^2} \int_{B(z,r)} \partial_i u(x,y) \, dx \, dy, $$
which follows from $\partial_i u$ itself being harmonic. By the divergence theorem, then this is equal to
$$ \frac{1}{\pi r^2} \int_{\partial B(z,r)} u \nu^i \, dS, $$
where $\nu^i$ is the $i$-th component of the unit vector normal to the surface $\partial B(z,r)$. Thus
$$ |\partial_i u(z)| \le \frac{1}{\pi r^2} \int_{\partial B(z,r)} |u| \, dS = \frac{2}{r} \|u\|_{L^{\infty}(\partial B(z,r))}$$
This allows you to control the pointwise convergence of the partial derivatives of $u_n$ by the uniform convergence of $u_n$. However, this in turn allows you to control the uniform convergence of the derivatives on a compact set, since then you can just use larger balls.
The reason your book uses the expressions $\frac{\partial f}{\partial z}$ and $\frac{\partial f}{\partial \overline{z}}$ without proof is that they are not actual partial derivatives, they are only derivatives in a symbolic sense. They are formal expressions that can help in understanding and communicating the behavior of a complex function.
If we write $x = \frac12(z + \overline{z})$ and $y = \frac{1}{2i}(z - \overline{z})$, thus expressing $x$ and $y$ as functions of $z$ and $\overline{z}$, and think of $f(z) = f(x,y)$ as a function of two real variables, then by supposing the rules of calculus apply we obtain
$$
\frac{\partial f}{\partial z} = \frac 12 \left(\frac{\partial f}{\partial x} - i \frac{\partial f}{\partial y} \right)
$$
and
$$
\frac{\partial f}{\partial \overline{z}} = \frac 12 \left(\frac{\partial f}{\partial x} + i \frac{\partial f}{\partial y} \right).
$$
If $\frac{\partial f}{\partial \overline{z}} = 0$ for a particular function, then this is equivalent with saying that $f$ satisfies the Cauchy-Riemann equations. So in a formal sense we could say that an analytic function is a function of $z$ alone and independent of $\overline{z}$, thinking of them independent variables. (Ahlfors)
Best Answer
$$\frac{d f}{d z}(z_0)=\frac{\partial u}{\partial x}(z_0) + i\frac{\partial v}{\partial x}(z_0)$$.