The Wasserstein space, which is the space of probability measures endowed with the so-called (quadratic) Wasserstein distance coming from optimal transport, can formally be seen as a Riemannian manifold of infinite dimension. We propose, through a variational approach, a definition of harmonic mappings defined over a domain of R^n and valued in the Wasserstein space. As the latter has nonnegative curvature, we cannot rely on the theory of Koorevaar, Schoen and Jost about harmonic mappings valued in metric spaces and we use arguments based on optimal transport instead. We manage to recover a fairly satisfying theory which captures some key features of harmonicity.