Composition of a harmonic function with a holomorphic function is still harmonic
The answer is yes, and we only need $g$ to be holomorphic. One can prove this by directly computing the Laplacian of $f\circ g$ using the Chain Rule. I'd rather use $z$ and $\bar z$ than $x$ and $y$ for this purpose.
$$\Delta(f\circ g)=\frac{1}{4}(f\circ g)_{z\bar z} = \frac{1}{4}[(f_z\circ g) g']_{\bar z} = \frac{1}{4}(f_{z\bar z}\circ g) \overline{g'} g'= [(\Delta f)\circ g]|g'|^2=0$$
Let $\phi(u,v)$ be harmonic in $D$. Let $w=f(z)=u(x,y)+iv(x,y)$ be analytic in $D_0$ defining a mapping $D_0\to D$. We have $$\phi_x=\phi_u u_x+\phi_v v_x$$ $$\phi_y=\phi_u u_y+\phi_v v_y$$ $$\phi_{xx}=\phi_{uu}(u_x)^2+\phi_{uv} u_x v_x +\phi_u u_{xx} +\phi_{vv} (v_x)^2+\phi_{vu} v_x u_x +\phi_v v_{xx}$$ $$\phi_{yy}=\phi_{uu}(u_y)^2+\phi_{uv} u_y v_y +\phi_u u_{yy} +\phi_{vv} (v_y)^2+\phi_{vu} v_y u_y +\phi_v v_{yy}$$ $$\phi_{xx}+\phi_{yy}=[(u_x)^2+(v_x)^2][\phi_{uu}+\phi_{vv}]$$ because $u_{xx}+v_{yy}=0$, $v_{xx}+v_{yy}=0$, $u_xv_x=-u_yv_y$. Hence $\phi_{uu}+\phi_{vv}=0$ implies $\phi_{xx}+\phi_{yy}=0$. We conclude that $\phi(x,y)$ is a harmonic funtion.
The term "harmonic" will mean "real and harmonic" below. I will use the following:
1) Compositions of holomorphic functions are holomorphic.
2) The real and imaginary parts of a holomorphic function are harmonic.
3) If $D\subset \mathbb C$ is an open disc, and $u$ is harmonic in $D,$ then there exists $v$ harmonic in $D$ such that $u+iv$ is holomorphic in $D.$
1) is just the chain rule. 2) follows from the Cauchy-Riemann equations. If you are unfamiliar with 3), I can add a proof later. I'll assume it for now.
Thm: Suppose $\Omega_1, \Omega_2 \subset \mathbb C$ are open. Assume $u$ is harmonic on $\Omega_2$ and $f: \Omega_1 \to \Omega_2$ is holomorphic. Then $u\circ f$ is harmonic on $\Omega_1.$
Proof: Let $a \in \Omega_1.$ Because harmonicity is a local property, it suffices to show $u\circ f$ is harmonic in a neighborhood of $a.$
Choose an open disc $D$ centered at $f(a)$ contained in $\Omega_2.$ Then by 3) above, there exists a $v$ harmonic in $D$ such that $u+iv$ is holomorphic in $D.$ Let $\omega = f^{-1}(D).$ Then $\omega $ is an open neighborhood of $a$ contained in $\Omega_1,$ and $f:\omega \to D.$ By 1), $(u+iv)\circ f$ is holomorphic in $\omega.$ By 2), the real part of $(u+iv)\circ f$ is harmonic in $\omega.$ But this real part is precisely $u\circ f,$ and we're done.
Your question can be interpreted in the greater context of "maps preserving harmonic functions".
Definition Let $(M,g)$ and $(N,h)$ be Riemannian manifolds. A mapping $\Phi:M\to N$ is said to be a harmonic morphism if whenever $u:N\to\mathbb{R}$ is a harmonic function (solving $\triangle_h u = 0$ where $\triangle_h$ is the Laplace-Beltrami operator for the Riemannian metric $h$) the composition $u\circ \Phi$ is a harmonic function on $M$.
Theorem A mapping is a harmonic morphism if and only if it is a harmonic map which is weakly horizontally conformal.
(Don't worry too much about the undefined terms in the above theorem.)
Corollary If $M$ and $N$ have the same number of dimensions, then
- If dimension is 2, $\Phi$ is a harmonic morphism if and only if $\Phi$ is conformal.
- If the dimension is bigger than 2, $\Phi$ is a harmonic morphism if and only if $\Phi$ is a conformal map with a constant coefficient of conformality.
For reference, see this paper of Bent Fuglede's.