Indeed the axes are completely arbitrary. The problem is spherically symmetric (i.e. the hamiltonian commutes with generators of rotation, which are the components of angular momentum).
However, we have chosen axes to describe the system and all our results are given with respect to these axes. Further, we arbitrarily chose a full set of commuting operators, whose eigenvalues uniquely specify any eigenvector of the hamiltonian. These operators are usually taken to be $H, L^2, L_z$ with their eigenvalues labled by $n$, $l$, $m$ respectively. This makes the $z$-axis special to our description of the system, while this of course it is still just an artifact of our choice of coordinates.
When you apply electric or magnetic fields to the atom the system is no longer spherically symmetric (since the fields come from some direction) and then different dumbbell shapes react differently to fields coming from different directions and having different polarization (this should be intuitively obvious). Hence my conclusion would be: The non-spherical shapes are what you get when you find simultanous eigenvectors of above operators, the choice of which is not spherically symmetric. In a problem, which is no longer spherically symmetric the shapes gain physical significance because you can actually excite electrons to those precise orbitals with light coming from appropriate directions and having appropriate polarization.