Is it a coincidence or a convention that the spherical harmonics $Y_{\ell m}(\theta, \phi)$ are eigenstates of $L_z$?

The $z$-component of the angular momentum $L_z$ is an eigenstate of the spherical harmonics $Y_{\ell m}(\theta, \phi)$. The spherical harmonics $Y_{\ell m}(\theta, \phi)$ are also the eigenstates of the total angular momentum operator $\mathbf{L}^2$.
This is well known in quantum mechanics, since $[\mathbf{L}^2, L_z] = 0$, the good quantum numbers are $\ell$ and $m$.

Would it be possible to find another solution analogous to the spherical harmonics $Y_{\ell m}(\theta, \phi)$ such that $[\mathbf{L}^2, L_{x \text{ or } y}] = 0$?

Is the $z$ direction just chosen a posteriori, or is there a meaning as to why only its component shares the spherical harmonics as eigenstates?

• Just take the spherical harmonics but with the $\theta$ and $\phi$ measured with respect to your chosen axis. There are formulas that relate that with the spherical harmonics with respect to the z axis. – Javier May 7 '17 at 22:32

There is no global $z$ axis that permeates the universe. As such, since the choice of $z$ axis is necessarily arbitrary for every specific situation, there is nothing special about the choice of $L_z$ as the operator whose eigenstates we choose as a basis; for every given axis $\hat n$, there exists a common eigenbasis of $L^2$ and $L_{\hat{n}}$.
As far as spherical harmonics go, we normally choose to define them as the eigenstates of $L_z$ (as opposed to some other axis) because the spherical coordinate system bestows a special significance on it (we measure the polar angle $\theta$ from the $z$ axis, and the azimuthal angle $\phi$ around the $z$ axis), but that's about it. You could equally well build a spherical coordinate system around the $x$ axis, and the $Y_{lm}$ evaluated on those coordinates will be eigenfunctions of $L_x$.
Just take the spherical harmonics but with the $\theta$ and $\phi$ measured with respect to your chosen axis. There are formulas that relate that with the spherical harmonics with respect to the $z$ axis,
it's actually a lot simpler than that. There are indeed ways to relate the $z$-polar coordinate angles (call them $\theta_z$ and $\phi_z$) with the $x$-polar angles $\theta_x$ and $\phi_x$, but they are rather complicated and quite nonlinear. However, you don't actually need them: the reason for this is that if you evaluate the spherical harmonics on the unit sphere, so that $x=\sin(\theta)\cos(\phi)$, $y=\sin(\theta)\sin(\phi)$ and $z=\cos(\theta)$, then the harmonics perform a magical transformation from ugly trigonometric functions into simple homogeneous polynomials in $x$, $y$ and $z$ of degree $l$, so, for instance, $$Y_{4,2}(\hat{\mathbf r}) \propto (x+iy)^2\left( x^2 + y^2 - 6z^2\right), \tag{assuming r^2=1}$$ where the first factor is always of the form $(x+iy)^m$ and the rest is real-valued and a function of $x^2 + y^2$ and $z^2$.
If you want to get that harmonic around the $x$ axis, then you can just do the rotation on the cartesian components (which is easy, even for any arbitrary axis, because the rotation is a linear transformation, under which the homogeneous polynomials are closed) and then sub in the polar coordinates as you require.