Can we characterize isotropic states via representation theory?

Question

The only bipartite states $\rho$ invariant under any $U\otimes \bar U$ transformation are isotropic states. More precisely, as discussed e.g. in quant-ph/0109124 in Section 4.4, a state $\rho$ is invariant under any such transformation iff it is a mixture of maximally entangled and maximally mixed states: $$\rho = p \mathbb{P}_+ + (1-p) \frac{I_N\otimes I_N}{N^2},$$ for some $p\in[0,1]$, where $N$ is the dimension of the underlying spaces (assumed to be equal), and $\mathbb{P}_+\equiv |+\rangle\!\langle+|$ projection onto the maximally entangled state $|+\rangle\equiv\frac{1}{\sqrt N}\sum_{k=1}^N |k,k\rangle$.

This was shown by some of the same authors in quant-ph/9708015, though if I understand correctly the main ideas were already discussed before in (Werner 1989). A rough sketch of the proof as presented in quant-ph/9708015 is as follows:

1. Let $A$ be $U\otimes\bar U$ invariant (for all unitaries $U$), and consider its matrix elements on the standard canonical (product) basis. By direct examination of what the invariant condition entails for unitaries $U$ that only act nontrivially on a single basis element, we observe that the only possible nonzero matrix elements are those of the form $A_{nn,mm}$ and $A_{nm,nm}$ for some (possibly equal) indices $n,m$.
2. Investigating the consequences of the condition when $U$ is a two-element permutation, we further find that invariance implies the form $A=bB+cC+dD$ for some $b,c,d\in\mathbb{R}$, where $B\equiv\sum_{m\neq n}|mn\rangle \!\langle mn|$, $C\equiv\sum_{m\neq n}|mm\rangle\!\langle nn|$, and $D\equiv\sum_{m}|mm\rangle\!\langle mm|$.
3. Considering now unitaries of the form $U=\tilde U_2\oplus I_{N-2}$ for some two-dimensional rotation $\tilde U_2$, we see that $D$ is not actually invariant, and thus $d=0$.
4. Imposing the condition $\operatorname{tr}(A)=1$, using $\operatorname{tr}(B)=N(N-1)$ and $\operatorname{tr}(C)=0$, we find that $b=1/[N(N-1)]$.
5. Observe that we ended up with an isotropic state for some $p\in[0,1]$, using the explicit expressions $\mathbb{P}_+=\frac{1}{N}\sum_{m,n}|m,m\rangle\!\langle n,n|$ and $I=\sum_{m,n}|m,n\rangle\!\langle m,n|$.

I have no issues with the proof above, but it's quite laborious. Is there a more direct (if more abstract) approach to obtain this result? For example, by using representation theory to characterise the irreps of the $\mathbf U(N)$ representation $\eta(U)(\rho)\equiv (U\otimes\bar U)\rho (U\otimes\bar U)^\dagger$? My hope was that starting from the simple observation that isotropic states are invariant under such transformations, we have $$\rho = \int dU\, (U\otimes\bar U)\rho (U\otimes\bar U)^\dagger,$$ and being the latter the commutant of the representation, it can be split into a direct sum of projections of $\rho$ onto the different irreps.

Would this work? My main issue is that I'm not sure how to fully characterise the irreps of this particular representation of $\mathbf U(N)$ on the space $\mathbf{GL}(\mathbb{C}^N\otimes\mathbb{C}^N)$.