I don't have a reference, but this is easy to show directly: we have $ \| ST \|_{HS}^2 = \sum_{n = 1}^{\infty} \| STe_n \|^{2} \leq \| S\|^2 \sum_{n = 1}^{\infty} \| Te_{n} \|^{2} = \| S \|^2 \| T \|_{HS}^2,$ where $\| \cdot \|_{HS}$ is the Hilbert-Schmidt norm. This, by itself, however is not enough to prove that the set of Hilbert-Schmidt operators form an ideal - assuming you mean two-sided ideal, you also need to show, with $S$ and $T$ as above, that $TS$ is Hilbert-Schmidt, and that the sum of two Hilbert-Schmidt operators is again Hilbert-Schmidt (along with the perfectly obvious statement that the zero-operator is Hilbert-Schmidt). The latter is easy, \| S + S' \|_{HS}^2 = \sum_{n = 1}^{\infty} \| (S + S')e_n \|^2 \leq \sum_{n = 1}^{\infty} (\| Se_n \|^2 + \| S'e_n \|^2) = \|S\|_{HS}^2 + \| S'\|_{HS}^2 < \infty if $S$ and S' are Hilbert-Schmidt, but I expect the former to be more difficult.
Edit: The result follows from the fact that the set of Hilbert-Schmidt operators is $^*$-closed, i.e. if $T$ is H-S, then so is $T^*$. The proof is easy: if $S \in B(H)$ and $T$ is H-S, then $(TS)^* = S^* T^*$ is H-S by the above, and hence $TS$ is H-S also.
Another edit: I'll add that the fact that the class of H-S operators is $^*$-closed follows from the fact that \| S \|_{HS} = \| S^* \|_{HS}. This is a (small) part of exercise 17c) of section XVIII.9 in Real and Functional Analysis by Serge Lang, so it shouldn't be too difficult to show (however, at the moment, I have no idea how to do it; maybe a sign that I need to revisit some of this stuff).