Convex functions defined on metric spaces are pulled back to subharmonic ones by harmonic maps

If $u:\Omega \subset \mathbb{R}^d \to X$ is a harmonic map valued in a metric space $X$ and $E : X \to \mathbb{R}$ is a convex function, in the sense that it generates an EVI-gradient flow, we prove that the pullback $E \circle u : \Omega \to \mathbb{R}$ is subharmonic. This property was known in the smooth Riemannian manifold setting or with curvature restrictions on $X$, while we prove it here in full generality. In addition, we establish generalized maximum principles, in the sense that the $L^q$ norm of $E \circ u$ on $\partial\Omega$ controls the $L^p$ norm of $E \circ u$ in $\Omega$ for some well-chosen exponents $p \geq q$, including the case $p=q=+\infty$. In particular, our results apply when $E$ is a geodesically convex entropy over the Wasserstein space, and thus settle some conjectures of Brenier (Optimal transportation and applications (Martina Franca, 2001), volume 1813 of lecture notes in mathematics, Springer, Berlin, pp 91-121, 2003).