The isotropic harmonic oscillator

While studying Lie-algebras I read that the three-dimensional harmonic oscillator has an $SU(3)$ symmetry. I found this very unexpected; I thought it was ''obvious'' that the symmetry is only $SO(3)$.
The Hamiltonian of the three-dimensional harmonic oscillator is (setting $\hbar=m=\omega = 1$ for simplicity) $$\begin{equation*} H = \frac{1}{2} \vec{p}^2 + \frac{1}{2} r^2 \end{equation*}$$ with $r^2 = x^2 + y^2 + z^2$. The energy levels are of the form [1] $$\begin{equation*} E_n = n + 3/2 \end{equation*}$$ with $n \ge 0$ an integer. The degeneracy of the $n$-th level $d(n)$ is equal to the number of ways $n$ can be written as $n = n_1 + n_2 + n_3$ where $0 \le n_1, n_2, n_3$ are integers. For example,

• If $n=0$, $d(0)=1$ because one only has $0 = 0+0+0$ 
• If $n=1$, $d(1)=3$ because one has the three possibilities $$\begin{align*} 1 &= 1+0+0\\ 1 &= 0+1+0\\ 1 &= 0+0+1 \end{align*}$$
• The case $n=2$ is already more interesting. $d(2) = 6$ coming from $$\begin{align*} 2 &= 2+0+0\\ 2 &= 0+2+0\\ 2 &= 0+0+2\\ 2 &= 1+1+0\\ 2 &= 1+0+1\\ 2 &= 0+1+1 \end{align*}$$ Here one can calculate that the states transform under $SO(3)$ as the representation $\bf{5} + \bf{1}$. There is thus a degeneracy in the energy levels that is not explained by the $SO(3)$ symmetry.
• At the next level, there is even more degeneracy: $d(3) = 10$ coming from $$\begin{align*} 3 &= 3+0+0,\quad\text{3 permutations in total}\\ 3 &= 2+1+0,\quad\text{6 permutations in total}\\ 3 &= 1+1+1\end{align*}$$ Here one can calculate that the states transform under $SO(3)$ as the representation $\bf{7} + \bf{3}$.

Mathematically, it is very easy to construct the $SU(3)$ symmetry. Because the harmonic oscillator is three-dimensional, there are three independent creation- and annihilation operators with commutation relations $$\begin{align*} [a_i, a_j] &=0\\ [a_i^{\dagger}, a_j^{\dagger}] &=0\\ [a_i, a_j^{\dagger}] &=\delta_{ij} \end{align*}$$ In general, if $T_a$ are $3\times3$ matrices satisfying $$\begin{equation*} [T_a,T_b] = \sum_c f_{ab}^c T_c \end{equation*}$$ then the operators $Q_a = \sum_{ij} a_i^{\dagger} T_{a,ij} a_j$ satisfy the same algebra $$\begin{equation*} [Q_a,Q_b] = \sum_c f_{ab}^c Q_c \end{equation*}$$ If the matrices $T_a$ form the 3-dimensional representation of $SU(3)$, then the $Q_a$ form a representation of the $SU(3)$ algebra, and it is easy to see that they also commute with the Hamiltonian. The states for $n=2$ form the $SU(3)$ representation $\bf{6}$. The states for $n=3$ form the $SU(3)$ representation $\bf{10}$, and so on. All in all, the mathematics is easy, but I would not have guessed that there is this extra symmetry. All this is explained in chapter 14 of [2], and also in notes by Kirson [3].
References
[1] Griffiths, Introduction to Quantum Mechanics, exercise 4.39
[2] Georgi, Lie Algebras In Particle Physics
[3] Kirson, Isotropic harmonic oscillator