You don't need to, but you can.
Essentially, what you're doing is choosing between expressing your energy eigenstates as products of functions of spherical, cylindrical, or Cartesian coordinates. Any one of these options is equally good, while in other problems the calculations will get quite cumbersome in, e.g., Cartesian coordinates.
At the end of the day, you need to specify each point in space with three numbers. Which separation of variables you use is just a matter of whether you prefer $(r, \theta, \phi)$, $(\rho, \phi, z)$, or $(x, y, z)$.
To be a bit more specific, if you use spherical coordinates, the solution will be written in terms of, e.g., spherical harmonics and some radial function. If you use Cartesian, no spherical harmonics will appear. Instead, you'll get products of Hermite polynomials in each direction.
For an example a bit more explicit, I'll notice that the solution in spherical coordinates is given on Wikipedia (e.g., on this article). In there it is mentioned that the energy levels are given by $$E_n = \hbar \omega \left(n + \frac{3}{2}\right),$$ where $n$ is a non-negative integer and $\frac{1}{2} \mu \omega^2 = a$ ($\mu$ is the mass of the particle and I'm changing the choice of constants to that used in Wikipedia so the notation gets closer to the standard when dealing with quantum harmonic oscillators). Wikipedia also mentions that the degeneracy at energy level $n$ (i.e., the dimension of the eigenspace with energy $E_n$) is given by $$\frac{(n+1)(n+2)}{2}.$$
Let us reproduce this result in Cartesian coordinates. I'll just sketch the calculations, but they are done in more detail in Sec. 2.5 of Weinberg's Lectures on Quantum Mechanics. In Cartesian, the Hamiltonian can be written as \begin{align} H &= - \frac{\hbar^2}{2 \mu} \nabla^2 + \frac{1}{2}\mu \omega^2 r^2, \\ &= \left(- \frac{\hbar^2}{2 \mu} \frac{\partial^2}{\partial x^2} + \frac{1}{2}\mu \omega^2 x^2\right) + \left(- \frac{\hbar^2}{2 \mu} \frac{\partial^2}{\partial y^2} + \frac{1}{2}\mu \omega^2 y^2\right) + \left(- \frac{\hbar^2}{2 \mu} \frac{\partial^2}{\partial z^2} + \frac{1}{2}\mu \omega^2 z^2\right), \end{align} which is just three times the Hamiltonian to a one-dimensional QHO. Hence, the system can be treated as three independent QHOs, and hence the allowed energy levels are
\begin{align} E_{qrs} &= \hbar \omega \left(q + \frac{1}{2}\right) + \hbar \omega \left(r + \frac{1}{2}\right) + \hbar \omega \left(s + \frac{1}{2}\right), \\ &= \hbar \omega \left(q + r + s + \frac{3}{2}\right), \end{align} where the non-negative integers $q$, $r$, and $s$ are the quantum numbers for each of the independent 1D-QHOs. Notice then that the allowed energies are $$E_n = \hbar \omega \left(n + \frac{3}{2}\right),$$ just as we found in the spherical separation of coordinates. However, is the degeneracy right? Fix $n$. How many combinations of $q$, $r$, and $s$ lead to the name $n$? Following Weinberg, we notice we have $$q + r + s = n.$$ Since $n$ is fixed, this means that once we choose values for $q$ and $r$, $s$ will be given by $s = n - q - r$. The degeneracy then must be given by \begin{align} \sum_{q = 0}^{n} \sum_{r = 0}^{n - q} 1 &= \sum_{q = 0}^{n} (n - q + 1), \\ &= (n + 1)^2 - \frac{n (n+1)}{2}, \\ &= \frac{(n+1)(n+2)}{2}, \end{align} which is exactly the same result.
To show the different separations of variables only lead to different choices of basis will likely take a bit more calculation, but I hope this helps to illustrate the result.