Is there a theory of the type of maps between domains that preserve harmonic functions? For instance, in the 2-dimensional case, we know that conformal maps (or even just holomorphic ones) are such maps. But what about higher dimensions, and what about other maps in this class of harmonic-function preserving maps?
In 2-d, We know for instance that any function that is conjugate analytic also preserves harmonic functions, and therefore as is any analytic function of $\overline{z}$. But are there others?
Such maps are called harmonic morphisms. I recommend the book Harmonic Morphisms Between Riemannian Manifolds by Baird and Wood.
To your specific question: a harmonic morphism $f$ from $U\subset \mathbb C$ to $\mathbb C$ is either holomorphic or anti-holomorphic. Indeed, $f$ is a harmonic morphism if $$h_{z\bar z}=0 \implies (h\circ f)_{z \bar z}=0 \tag{1}$$ We can allow complex-valued $h$ here, since $h$ enters (1) linearly anyway. Using $h(z)=z$ yields $f_{z\bar z}=0$. Using $h(z) = z^2$ yields $$ (f^2)_{z \bar z} = 2 f f_{z \bar z} +2 f_z f_{\bar z} = 2f_zf_{\bar z} $$ Since $f_z$ is holomorphic and $f_{\bar z}$ is antiholomorphic, both of them are either identically zero or have a discrete zero set. Since $f_zf_{\bar z}\equiv 0$, it follows that one of the factors is identically zero.