The worst you can do, is the 'democratic' coupling $\pi := \sum_{i,j} a_ib_i\delta_{x_i}\otimes \delta_{y_j}$, which gives you the upper bound $$M= ( \sum_{i,j} a_ib_j |x_i-y_j|^p )^{\frac{1}{p}}.$$

Okay, I was making first sure to identify your question. I did not report a proof in the answer, but links to read this fact. It is all in the references I gave, e.g. (a)-(b) in Theorem 4.16 of arxiv.org/abs/1106.2090.

For $p=\infty$, there is a 'weak maximum principle' to be shown, namely $$ f \le C, \quad a.e. \quad \Rightarrow \quad h_t(f) \le C, \quad a.e.,\forall t>0.$$ Again, this can be shown to be true to cover the case of $L^\infty$

One can prove that, for $u \colon \mathbb{R}\rightarrow [0,\infty]$ convex and l.s.c. with $u(0)=0$, the mapping $t\mapsto \int u(h_t(f))\, dVol_g$ is nonincreasing. Apply this to $u(\cdot)=|\cdot|^p$, for any $p \in [1,\infty)$ to get $$\|h_t(f)\|_{L^p}\le \|f\|_{L^p}, \qquad \forall f \in L^2\cap L^p, t\ge 0.$$ The case $p=\infty$ is to be treated differently. Is this what you asked?

If you do not want to assume the existence of a heat kernel, you may consider the heat semigroup as a gradient flow associated with the Dirichlet energy. In this framework, the extension to L^p spaces and contractivity is a standard result in the theory and holds true even in the setting of the so-called metric measure spaces, containing the class of Riemannian Manifolds. Would be this setting good to cite for you?

Absolutely right: either you fix $\mu_{0,n},\mu_{1,n}$ and $t\mapsto \mu_{t,n}$ as hypothesis of your problem, and and try to exhibit $\mu_t$, or you impose some additional properties on the base space $X$ (e.g. Geodesic Space) not to deal with existence issues. It depends on which statement you seek.

I think is very hard to fix $\mu_t$ and produce, afterwards, approximating marginals $(\mu_{0,n}),(\mu_{1,n})$. If you stick to the other way, as in my answer, you can try to argue essentially by tightness to get $\mu_t$ from $\mu_{t,n}$. Finally, you just need to show that the limit $\mu_t$ is a $W_2$ geodesic. The first link I gave you follows this path. With this approach, is very hard to control whether you are converging to your fixed a priori Wasserstein geodesic, or to another.

This is a comment rather than an answer, but I could not post it as a comment. Anyway, something useful in this direction can be found in Lemma 4.4 arxiv.org/pdf/1609.00782.pdf which, combined with Proposition 4.8 of arxiv.org/pdf/1311.4907.pdf gives you $W_2$ close $\mu_{t,n}$'s.