In this article, we introduce a new approach towards the statistical learning problem $\operatorname{argmin}_{ρ(θ) \in \mathcal P_θ} W_{Q}^2 (ρ_{\star},ρ(θ))$ to approximate a target quantum state $ρ_{\star}$ by a set of parametrized quantum states $ρ(θ)$ in a quantum $L^2$-Wasserstein metric. We solve this estimation problem by considering Wasserstein natural gradient flows for density operators on finite-dimensional $C^*$ algebras. For continuous parametric models of density operators, we pull back the quantum Wasserstein metric such that the parameter space becomes a Riemannian manifold with quantum Wasserstein information matrix. Using a quantum analogue of the Benamou-Brenier formula, we derive a natural gradient flow on the parameter space. We also discuss certain continuous-variable quantum states by studying the transport of the associated Wigner probability distributions.