Let ${\cal X }=XX^{\prime}$ be a random matrix associated with a centered $r$-column centered Gaussian vector $X$ with a covariance matrix $P$. In this article we compute expectations of matrix-products of the form $\prod_{1\leq i\leq n}({\cal X } P^{v_i})$ for any $n\geq 1$ and any multi-index parameters $v_i\in\mathbb{N}$. We derive closed form formulae and a simple sequential algorithm to compute these matrices w.r.t. the parameter $n$. The second part of the article is dedicated to a non commutative binomial formula for the central matrix-moments $\mathbb{E}\left(\left[{\cal X }-P\right]^n\right)$. The matrix product moments discussed in this study are expressed in terms of polynomial formulae w.r.t. the powers of the covariance matrix, with coefficients depending on the trace of these matrices. We also derive a series of estimates w.r.t. the Loewner order on quadratic forms. For instance we shall prove the rather crude estimate $\mathbb{E}\left(\left[{\cal X }-P\right]^n\right)\leq \mathbb{E}\left({\cal X }^n-P^n\right)$, for any $n\geq 1$