4D chemistry is funny. Yes, that was a great introductory paragraph.
So, if you try to solve the Schrodinger equation in 4D, then the result is that the radial part (how the wavefunction changes with distance) is a doofus, but the angular part is well behaved. So to analyze orbitals, the best way is to analyze the radial part.
In 3D, the radial part is a spherical harmonic. This means it is an eigenfunction of the squared angular momentum operator. In 4D, the radial part will be a hyperspherical harmonic. This means it is an eigenfunction of the squared angular momentum operator. (who could have guessed?)
We can restrict the radial part some more. In 3D, we can measure both the angular momentum about the z-axis and the squared angular momentum simultaneously. So we can demand that the spherical harmonics are eigenfunctions of squared angular momentum and angular momentum about the z-axis. It turns out that in 4D, we can force the hyperspherical harmonics to be eigenfunctions of squared angular momentum, angular momentum in the x-y plane, and angular momentum in the z-w plane.
The SO(4) algebra
The angular momentum operator in the $x_i$-$x_j$ plane is $L_{ij}=i\left(x_i\frac{\partial}{\partial x_j}-x_j\frac{\partial}{\partial x_i}\right)$. Or maybe the negative of that. The sign isn’t important, for reasons that may or may not become clear later.
We can consider the operators $[L_{ij},L_{kl}]=L_{ij}L_{kl}-L_{kl}L_{ij}$, known as commutators. These are important for a variety of reasons. For example, if $[L_{ij},L_{kl}]=0$, then they can be simultaneously measured, and otherwise, Heisenberg’s uncertainty principle applies to them.
The commutators also allow us to analyze the hyperspherical harmonics. If we let $Y^\ell_{m,n}$ be the hyperspherical harmonic with squared angular momentum $\ell(\ell+2)$, angular momentum in the x-y plane $m$, and angular momentum in z-w plane $n$, then we can apply so-called ladder operators to derive other hyperspherical harmonics.
We know that $L_{12}Y^\ell_{m,n}=m$. Suppose that $J$ was a ladder operator. Then $L_{12}JY^\ell_{m,n}=JL_{12}Y^\ell_{m,n}+[L_{12},J]Y^\ell_{m,n}=(mJ+[L_{12},J])Y^\ell_{m,n}$. Note that if $[L_{12},J]=\pm J$, then $L_{12}JY^\ell_{m,n}=(m\pm 1)JY^\ell_{m,n}$. Similarly, if $[L_{34},J]=\pm J$, then $L_{34}JY^\ell_{m,n}=(n\pm 1)JY^\ell_{m,n}$. So if both $[L_{12},J]=\pm J$ and $[L_{34},J]=\pm J$, then $JY^\ell_{m,n}$ must be a multiple of $Y^\ell_{m\pm 1,n\pm 1}$.
By messing around, we obtain the following four operators:
\begin{gather*}
J_{++}=L_{13}-iL_{23}-iL_{14}-L_{24}\\
J_{+-}=L_{13}-iL_{23}+iL_{14}+L_{24}\\
J_{-+}=L_{13}+iL_{23}-iL_{14}+L_{24}\\
J_{- -}=L_{13}+iL_{23}+iL_{14}-L_{24}
\end{gather*}
These are named so that, e.g. $J_{+-}$ increases $m$ by 1 and decreases $n$ by 1.
Stuff with orbitals
With these operators, note that $Y^\ell_{0,0}$ and $Y^\ell_{1,0}$ are independent: no application of $J$s can send $Y^\ell_{0,0}$ to $Y^\ell_{1,0}$. There may also be the orbitals obtained from $Y^\ell_{\frac 12,\frac 12}$ and $Y^\ell_{\frac 12,-\frac 12}$. These would seem to create two (or four) classes of orbitals.
However, it turns there is only one class of orbital. This is generated from $Y^\ell_{-\ell,0}=(x+iy)^\ell$ (in some arbitrary normalization I just specified). Applying the $J_{\pm\pm}$ operators generates a family of $\ell^2$ hyperspherical harmonics. It can be shown that for angular momentum $\ell$, there are indeed $\ell^2$ hyperspherical harmonics, so these functions are all the hyperspherical harmonics.
Therefore, the shells are labeled by principal quantum number $n$ and angular momentum $\ell$. Each shell holds $4\ell^2$ electrons (4 because in 4D, electrons have 4 spins). It’s reasonable to suppose they fill in order of increasing $(n+\ell,n)$, as there isn’t any way to determine the actual filling order.