Where do symmetries in atomic orbitals come from? It is well established that: 
'In quantum mechanics, the behavior of an electron in an atom is described by an orbital, which is a probability distribution rather than an orbit. 
There are also many graphs describing this fact:
http://en.wikipedia.org/wiki/Electron:
(hydrogen atomic orbital - one electron)
In the figure, the shading indicates the relative probability to "find" the electron, having the energy corresponding to the given quantum numbers, at that point.
My question is: How do these symmetries shown in the above article occur?
                What about the 'preferable' axis of symmetries? Why these?
 A: 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. 
A: 
How do these symmetries shown in the above article occur? What about the 'preferable' axis of symmetries? Why these?

For atoms subject to no net external electric of magnetic fields the orientation of the axes is arbitrary. This shows up clearly in the math because adding up all the spherical harmonic contributing to a single shell (1s, 2s, 2p, 3s, 3p, 3d, ...) gives no angular dependence. It doesn't show clearly in the visualization because those plots employ an arbitrary cut-off in generating the display. So, short answer, the lobes of the orbitals point along the coordinate axes purely for convenience: there is no physics content to that feature of the rendering.
The fact that there are a non-negative integer number of radial or angular nodes arises from the boundary conditions on the wave-function: just like the vibrations of a guitar string only those modes that 'fit' in the space exist as time-independent solutions.
In the case that there are external electromagnetic fields, then those fields do two things:


*

*They change the shape of the time-independent solutions

*The enforce a choice of orientation on the new solutions.

A: The short answer is the electrons need to obey rules when they are bound to the atom.
I orginally had a picture of orbitals here, but please follow dmckee's comments, in which the pitfalls of overreliance on pictures, rather than math visualisation, is correctly pointed out. If you can work at the math, it will eventually reward your effort. 


*

*They occupy the lowest energy levels first and then fill outwards with increasing energy.

*They follow the Pauli Exclusion Principle, so two electrons can only occupy the same orbital if they have opposite spins.

*They are mutually repulsive, because of their negative charge, so they may not be symmetrical in space, but they will be as far apart as they can get from one another. They are symmetric in the sense that seen in a mirror, they would appear the same, although not spherically symmetric as this may bring them closer together and break the rules that govern their "location" 

*Each higher energy level may allow a lot more electrons to occupy it, in various shapes, but once an energy level is full, any new arrival of anelectron must occupy a discrete orbital, and all of these factors together takes us further away from the orginal Bohr solar system model.
A: 
Where do symmetries in atomic orbitals come from?

It is always a good idea to search the historical development. I recommend you to read the Wikipedia article about Spherical harmonics. Spherical harmonics are special functions defined on the surface of a sphere. They were first investigated in connection with the Newtonian potential of Newton's law of universal gravitation in three dimensions by Pierre-Simon de Laplace in 1782.

The prevalence of spherical harmonics already in physics set the stage for their later importance in the 20th century birth of quantum mechanics. The spherical harmonics are eigenfunctions of the square of the orbital angular momentum operator ... and therefore they represent the different quantized configurations of atomic orbitals.


Visual representations of the first few real spherical harmonics. Blue portions represent regions where the function is positive, and yellow portions represent where it is negative (from Wikipedia https://en.wikipedia.org/wiki/Spherical_harmonics#/media/File:Spherical_Harmonics.png)

What about the 'preferable' axis of symmetries? Why these?

As you can see in the following image from the same Wikipedia article the orientation of the spherical harmonics to the axis X,Y and Z of the Cartesian coordinates are not the only possible orientations:

From Wikipedia: Sperical harmonics
As you can see from this figures the symmetrical figures due to Cartesian coordinates are the usual way of thinking but not the only possible.
