is there a way to prove the Mean-value formulas using complex analysis? 
THEOREM : $U$ is an open subset of $\mathbb{R^n}$ and suppose $u \in C^2(U) $ is harmonic within $U$, then :
  $$u(x)= \def\avint{\mathop{\,\rlap{-}\!\!\int}\nolimits} \avint_{B(x,r)}u\,dy = \def\avint{\mathop{\,\rlap{-}\!\!\int}\nolimits} \avint_{\partial B(x,r)}u\,dS$$
  for each ball $B(x,r) \subset U$

there's a nice proof in Evan's PDE book which only involves real analysis. since harmonic functions can be implicitly identified as the real or imagniary part of holomorphic functions, how does one prove or atleast reformulate this theorem using complex analysis notions ?
any comment, references will be greatly appreciated.  
 A: For $n=2$ one can proceed as follows (this is essentially an elaboration
of above comments):
If $u$ is harmonic in $ U \subset \Bbb C$ and $B(z_0, R) \subset U$ then
$u = \operatorname{Re} f$ for some holomorphic function $f$ in $B(z_0, R)$.
This is true in any simply-connected domain in $\Bbb C$, $f$ can for
example be chosen as an anti-derivative of the (holomorphic) function
$u_x - i u_y$.
The Cauchy integral formula then states that for $0 < r < R$ and the path
$\gamma(t) = z_0 + re^{it}$, $0 \le t \le 2 \pi$,
$$
 f(z_0) = \frac{1}{2 \pi i} \int_\gamma \frac{f(\zeta)}{\zeta - z_0} \, d\zeta = \frac{1}{2 \pi} \int_0^{2 \pi} f(z_0 + re^{it}) \, dt \, .
$$
Taking real parts gives the mean-value formula for $u$:
$$
 u(z_0) =  \frac{1}{2 \pi} \int_0^{2 \pi} u(z_0 + re^{it}) \, dt \, .
$$
Finally, the “area form” of the mean-value formula follows by integration
over the radius:
$$
 \frac{1}{\pi r^2} \int_{B(z_0, r)} u(x, y)\, dx dy =
 \frac{1}{\pi r^2} \int_0^r \int_0^{2 \pi}  u(z_0 + \rho e^{it})
 \rho \, dt \, d\rho = \frac{1}{\pi r^2} \int_0^r 2 \pi u(z_0) \rho \, d\rho = u(z_0) \, .
$$
