The hydrogen atom is spherically symmetric, so for any solution of the Schrödinger equation for the hydrogen atom, any rotation of that solution must also be a solution. If you do the math on how to rotate a solution, it turns out that the solutions with a particular energy $E_n$ fall into groups labeled by an integer $l < n$. The integer $l$ is physical: $\hbar^2 l(l+1)$ is the magnitude squared of the angular momentum. Within each group, rotating the solution gives you a new solution in the same group. These two facts are of course connected: a rotation can't change the length of a vector.
One can show that each group contains $2l+1$ independent solutions, in that any solution $|n,l\rangle$ where the energy is $E_n$ and the angular momentum $\hbar^2 l(l+1)$ can be written as a sum $$|n,l\rangle = \sum_{m=-l}^l c_m |n,l,m\rangle$$ (I apologize for the somewhat poor notation.)
This decomposition is based on choosing a particular axis, and taking each state to depend on the angle $\varphi$ around this axis as $e^{im\varphi}$. The appearance of axes of symmetry in these plots is due to this choice of axis and particular decomposition. With another choice of axis, which is the same as a rotation, the states will be mixed.
The bottom line is that it's not each solution -- wavefunction -- that needs to be spherically symmetric, but the total set of solutions.