Why are hydrogen energy levels degenerate in $\ell$ and $m$? Is there a good physical picture of why the energy levels in a hydrogen atom are independent of the angular momentum quantum number $\ell$ and $m$?
 A: The shortest and correct answer: this degeneracy is determined by the symmetry of the system.
The case of degeneracy in hydrogen atom is so-called "accidental degeneracy", when eigenfunctions belonging to different irreducible representations of the symmetry group of a Hamiltonian correspond to the same energy. This type of degeneracy can also occurs in larger systems, for instance, in molecules. This degeneracy can not be predicted only from the standart consideration of  Hamiltonian symmetry. The reason for this degeneracy is the existance of hidden symmetry in the system.
Mathematically it means that one can construct for the systems with hidden symmetry some conserved quantities, so-called "integral invariant", which should be included in consideration of symmetry propetries in addition to symmetry of Hamiltonian. And in principle, it is possible to solve Schroedinger equation in more tricky manner with inclusion of these "integral invariants" and obtain solution, for which these "accidental degeneracies" will be strictly included.
In the case of hydrogen atom the reason is the invarience of the system not only to the three-dimensional rotation group 0(3), but also to the four-dimensional rotation group 0(4) - the system have unexpected 
at first glance hidden symmetry.
A: The degeneracy of energy-levels can be traced to the fact that the hydrogen atom possesses an enhanced $SO(4)$ symmetry caused by (among other things) the conservation of the Laplace-Runge-Lenz vector operator, see e.g. this Phys.SE post and Ref. 1.
References:


*

*G. 't Hooft, Introduction to Lie Groups in Physics, lecture notes, chapter 9. The pdf file is available here.

A: Supersymmetry
Well, for fixed $l$, the degeneracy of $m$ is because of SO(3) symmetry, we are just seeing a full representation of this group.  
The big question is why all the radial hamiltonians $H_l$ for different angular momenta have the same spectrum except a discrete number of eigenvalues.
Note that particularly the tower-spectrum for $l$ and the tower for $l+1$ only differ in one eigenvalue, the lowest energy one. This is the typical setup one can read in Witten's Supersymmetric Quantum Mechanics: a pair of hamiltonians differing only in the vacuum eigenstate. So you should be able to build a supersymmetry generator Q such that $H_1=QQ^+$ is the radial hamiltonian for angular momentum $l$ and $H_2=Q^+Q$ is the radial hamiltonian for angular momentum $l+1$. 
SUSY QM is simpler than QFT QM; it does not contemplate Spin; the state and the superpartner are just two levels in QM hamiltonians. It is just a little bit more advanced, mathematically, that the factorisation method; still it allows for some topological arguments on susy breaking that generalize to the QFT version, this was the idea of Witten when defining it. 
Just now I am not sure if this the connection for every potential having $l,m$ degeneracy does exist, or only for Coulomb-Hydrogen case; to start with, it implies that the potential $V(r)$ must come from a superpotential, so surely it is not so trivial to do, not to classify all the families of radial potentials that allow to do this trick. But is is a twenty years old idea by now, so surely it is already done.
Ok, even there is an entry in the wikipedia. According it, the superpotential is
$$W = \frac{\sqrt{2m}}{h} \frac{\lambda}{2(l+1)} - \frac{h(l+1)}{r\sqrt{2m}}$$
So that the potentials
$$V_-=W^2-W'= -\lambda \frac{1}{r} + \frac{h^2 l (l+1)} {2m} \frac{1}{r^2}- \frac{\lambda^2 m}{2  h^2 (l+1)^2} $$
and
$$V_+=W^2+W'=-\lambda \frac{1}{r} + \frac{h^2 (l+1) (l+2)} {2m} \frac{1}{r^2} + \frac{\lambda^2 m}{2  h^2 (l+1)^2}$$
have the same spectrum except for the lowest energy eigenvalue of the first one, which is zero and can not appear in the second one (nice topological result).
The advantage of this explanation is that it can be extended to potentials without the full SO(4) symmetry and to more exotic cases where the pairing fails for other eigenvalues.
PS: It can be noticed that the superpotential for Coulomb problem is just a constant away from $W(r)=1/r$. An interesting point is that this superpotential can be calibrated to pair with the free particle: $V_+(r)=W^2+W'=0$; the superpotentials having this property generate the so-called "transparent potentials", with special properties in the phase-shift. They can be thought as generalizing the radial equation to symmetric spaces, with $1/r$ being the euclidean case.
connection with so(4) group representations (and Runge-Lenz vector?)
According last page of this lecture, the role of Runge-Lenz vector as a supercharge is analytically tricky. But at least we get some help from group theory, if we recall that so(4) ~ su(2) x su(2) and that so(3) ~ su(2). So for our purposes we could really write
$$so(4) \approx su(2) \oplus so(3)$$
The rotational part, so(3), gives us the $m$ degeneracy inside a representation of the group of rotations; this should exists for every central potential. The $su(2)$ part is the one whose ladder operator allows an interpretation as supersymmetry charge, where the degeneracy is not complete because of the difference in the lowest energy eigenstate; usually it disappears because $Q |\Omega>$ is not normalizable. (I am a bit fascinated that the susy ladder operator is related to an SU(2), because in susy QM this is not needed, or at least not explicit)
The supersymmetry generator Q (alt. $Q^+$) when applied to the eigenfunction of a hamiltonian produces the corresponding eigenfunction in the partner. This is the same role that the ladder generator used to produce states inside a representation of the symmetry group, but in this view it comes from the susy pairing: if
$$H_2 \psi = Q^+Q \psi = E\psi$$
then
$$H_1 (Q\psi) = (Q Q^+) (Q \psi) = Q H_2 \psi = E (Q\psi)$$
Note that the pairing fails if $Q\psi$ doesn't exist; this is the case for the vaccum potential, but I remember that J Casahorran did some study for other eigenstates beyond the vacuum (it is tricky because of Witten's results). 
