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

If $u : \Omega\subset \mathbb{R}^d \to {\rm X}$ is a harmonic map valued in a metric space ${\rm X}$ and ${\sf E} : {\rm X} \to \mathbb{R}$ is a convex function, in the sense that it generates an ${\rm EVI}_0$-gradient flow, we prove that the pullback ${\sf E} \circ u : \Omega \to \mathbb{R}$ is subharmonic. This property was known in the smooth Riemannian manifold setting or with curvature restrictions on ${\rm 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 ${\sf E} \circ u$ on $\partial \Omega$ controls the $L^p$ norm of ${\sf 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 ${\sf E}$ is a geodesically convex entropy over the Wasserstein space, and thus settle some conjectures of Y. Brenier, "Extended Monge-Kantorovich theory" in Optimal transportation and applications (Martina Franca, 2001), volume 1813 of Lecture Notes in Math., pages 91-121. Springer, Berlin, 2003.