A question about orthogonal functions and Laplace's equation

Question

Suppose we have a harmonic function in spherical coordinates, $$ \nabla^2 f(r,\theta,\phi) = 0 $$ By using the separation of variables method to solve the PDE, we will find $$ f(r,\theta,\phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^l \left( A_l^m r^l + B_l^m r^{-l-1} \right) Y_l^m(\theta,\phi) $$ and it so happens that the spherical harmonic $Y_l^m(\theta,\phi)$ functions also form a complete orthogonal basis for the sphere, so that $$ g(\theta,\phi) = \sum_{lm} g_l^m Y_l^m(\theta,\phi) $$ Now my question is suppose I want to change my coordinate system to some $\alpha(r,\theta,\phi),\beta(r,\theta,\phi),\gamma(r,\theta,\phi)$. If I then solve Laplace's equation in this new system, $$ \nabla^2 f(\alpha,\beta,\gamma) = 0 $$ will I also find some sort of functions similar to $Y_l^m(\theta,\phi)$ which form a complete orthogonal basis for some 2D subspace? Or is spherical coordinates special in the sense that this wouldn't happen in just any coordinate system?