Radial Fourier transforms provide a good, consistent perspective on most of the theory. The Fourier transform $\widehat{f}(t)$ of a function $f \colon \mathbb{R}^n \to \mathbb{R}$ is given by the integral of $f(x) e^{2\pi i \langle x,t \rangle} \, dx$ over $x \in \mathbb{R}^n$. If $f$ is a radial function (i.e., $f(x)$ depends only on $|x|$), then we can radial symmetrize everything and the exponential function averages out to a radial function. Specifically, we get $$ \widehat{f}(t) = 2\pi |t|^{-(n/2-1)} \int_0^\infty f(r) J_{n/2-1} (2 \pi r |t|) r^{n/2} \, dr. $$ The precise factors are a little annoying, but basically this just means $J_{n/2-1}$ is what you get when you radially symmetrize an exponential function in $n$ dimensions. It's easy to see that if you symmetrize $e^{2\pi i \langle x,t \rangle}$ by averaging over all $x$ on a sphere, then you get a radial function of $t$, and furthermore as you vary the radius of the sphere you just rescale the function. So the one function $J_{n/2-1}$ captures all of this, modulo scaling.
One consequence is that Bessel functions inherit the orthogonality of the exponential functions (i.e., the different scalings are orthogonal), so they also inherit all the consequences of orthogonality. For example, this is really where the differential equation comes from. There's a strong analogy between Bessel functions and orthogonal polynomials, where rescaling the Bessel function corresponds to varying the degree of the polynomial.
You also get certain qualitative results for free: for example, the product of two Bessel functions should be an integral of Bessel functions with positive coefficients, since this corresponds to saying the product of two radial, positive-definite functions remains positive definite. You can write down the coefficients explicitly, but sometimes all you need is nonnegativity, and in any case this point of view makes it easy to believe that there should be an explicit formula.
This is basically a low-brow version of the representation theory approach. Basically, ordinary Fourier analysis studies $L^2(\mathbb{R}^n)$ under the action of the translation group $\mathbb{R}^n$. If you look at the full group of isometries of $\mathbb{R}^n$ (including the orthogonal group), then it's just a little more elaborate, and the Bessel functions arise as zonal spherical functions. It's worthwhile working through this perspective, but in practice just thinking about radial Fourier analysis gives you most of the benefits with less machinery.

