what else is in $\prod_{j=1}^n(1+q^j)$?
Up to scaling, $\prod_{j=1}^n(1+q^j)$ is the character of the principal specialization of the spinor representation of $\mathfrak{so}(2n+1)$. This was first explicitly stated by J. W. B. Hughes, Lie algebraic proofs of some theorems on partitions, in Number Theory and Algebra, Academic Press, 1977, pp. 135--155.
In a $2n$ dimensional vector space equipped with a symplectic form, you can look at all the Lagrangian subspaces, and this forms a subvariety of the Grassmannian $G(n,2n)$ called the Lagrangian Grassmannian, $\Lambda(n)$. It is a smooth projective variety which inherits a schubert cell decomposition from the Grassmannian, and lots of other Schubert calculus goodness. It's cells are indexed by strict partitions with largest part at most $n$. The number of points of $\Lambda(n)$ over finite fields is given by $\prod_{i=1}^n(1+q^i)$, so unimodality of the coefficients can be obtained from the Hard Lefschetz theorem. (Not like this is a simpler way of looking at it, but perhaps a geometric point of view can be useful for certain applications, for example the classical q-binomial theorem follows from a certain stratification of $\Lambda(n)$.)
A reference for an explicit presentation of the cohomology ring using Schur Q-polynomials can be found in "Algebro—Geometric applications of Schur S- and Q-polynomials" by Piotr Pragacz, Topics in Invariant Theory ,M. P. Malliavin Ed., Springer Lecture Notes in Math., 1478, 1991.
Let $exp(z,q)=\sum_{k=0}^{\infty}z^k/[k]_q!$ be the usual $q-$analogue of the exponential function $e^x.$
The identity
$\sum_{k=0}^n{q^k}{\binom{n}{k}_{q^2}}=\prod_{j=1}^n(1+q^j)$
can be obtained by comparing coefficients in
$$exp(\frac{z}{1+q},q^2)exp(\frac{qz}{1+q},q^2)=exp(z,q),$$
which is a natural $q-$analogue of $e^\frac{x}{2} e^\frac{x}{2}=e^x.$
This identity occurs in Séminaire Lotharingien de Combinatoire, B05a (1981), but is perhaps older.