[proofplan]
We first check that the interpolated measures have finite $p$-moment, using the fact that every path in $\operatorname{Geo}(X)$ is controlled by its two endpoints. For two times $s,t$, the pushforward $(e_s,e_t)_\#\Pi$ gives an admissible transport plan between $\mu_s$ and $\mu_t$, and the constant-speed property of each path bounds its transport cost by $|t-s|^p$ times the endpoint cost. Since the endpoint plan is optimal, this gives an upper bound for every intermediate distance. The triangle inequality for $W_p$ along the chain $\mu_0,\mu_s,\mu_t,\mu_1$ forces all these upper bounds to be equalities.
[/proofplan]
[step:Verify that every interpolated measure has finite $p$-moment]
Fix a point $x_\ast\in X$. For each $t\in[0,1]$, the map $e_t:\operatorname{Geo}(X)\to X$ is Borel, so $\mu_t=(e_t)_\#\Pi$ is a Borel probability measure on $X$.
We prove that $\mu_t$ has finite $p$-moment. Since each $\gamma\in\operatorname{Geo}(X)$ is a constant-speed geodesic,
\begin{align*}
d(\gamma(0),\gamma(t))=t\,d(\gamma(0),\gamma(1)).
\end{align*}
By the triangle inequality in $X$,
\begin{align*}
d(x_\ast,\gamma(t))\le d(x_\ast,\gamma(0))+t\,d(\gamma(0),\gamma(1)).
\end{align*}
Again by the triangle inequality in $X$,
\begin{align*}
d(\gamma(0),\gamma(1))\le d(\gamma(0),x_\ast)+d(x_\ast,\gamma(1)).
\end{align*}
Therefore, since $0\le t\le 1$,
\begin{align*}
d(x_\ast,\gamma(t))\le 2d(x_\ast,\gamma(0))+d(x_\ast,\gamma(1)).
\end{align*}
Using the elementary inequality $(a+b+c)^p\le 3^{p-1}(a^p+b^p+c^p)$ for $a,b,c\ge0$, we get
\begin{align*}
d(x_\ast,\gamma(t))^p\le 3^{p-1}\left(2^p d(x_\ast,\gamma(0))^p+d(x_\ast,\gamma(1))^p\right).
\end{align*}
Integrating with respect to $\Pi$ and using $(e_0,e_1)_\#\Pi=\pi$, whose marginals are $\mu_0$ and $\mu_1$, gives
\begin{align*}
\int_X d(x_\ast,z)^p\,d\mu_t(z)\le 3^{p-1}\left(2^p\int_X d(x_\ast,x)^p\,d\mu_0(x)+\int_X d(x_\ast,y)^p\,d\mu_1(y)\right).
\end{align*}
The right-hand side is finite because $\mu_0,\mu_1\in\mathcal P_p(X)$. Hence $\mu_t\in\mathcal P_p(X)$.
[/step]
[step:Use the dynamical plan to bound the distance between two time slices]
Let $s,t\in[0,1]$ with $s\le t$. Define the endpoint-at-times coupling
\begin{align*}
\pi_{s,t}:=(e_s,e_t)_\#\Pi.
\end{align*}
Then $\pi_{s,t}\in\Pi(\mu_s,\mu_t)$ because its first marginal is $(e_s)_\#\Pi=\mu_s$ and its second marginal is $(e_t)_\#\Pi=\mu_t$.
By the definition of $W_p$ as the infimum of $p$th transport costs over couplings,
\begin{align*}
W_p(\mu_s,\mu_t)^p\le \int_{X\times X} d(x,y)^p\,d\pi_{s,t}(x,y).
\end{align*}
Using the pushforward identity for $\pi_{s,t}$, this becomes
\begin{align*}
W_p(\mu_s,\mu_t)^p\le \int_{\operatorname{Geo}(X)} d(\gamma(s),\gamma(t))^p\,d\Pi(\gamma).
\end{align*}
Since each $\gamma\in\operatorname{Geo}(X)$ is parametrized with constant speed on $[0,1]$,
\begin{align*}
d(\gamma(s),\gamma(t))=(t-s)d(\gamma(0),\gamma(1)).
\end{align*}
Therefore
\begin{align*}
W_p(\mu_s,\mu_t)^p\le (t-s)^p\int_{\operatorname{Geo}(X)} d(\gamma(0),\gamma(1))^p\,d\Pi(\gamma).
\end{align*}
Using $(e_0,e_1)_\#\Pi=\pi$, we identify the last integral with the endpoint transport cost:
\begin{align*}
\int_{\operatorname{Geo}(X)} d(\gamma(0),\gamma(1))^p\,d\Pi(\gamma)=\int_{X\times X} d(x,y)^p\,d\pi(x,y).
\end{align*}
Since $\pi$ is optimal,
\begin{align*}
\int_{X\times X} d(x,y)^p\,d\pi(x,y)=W_p(\mu_0,\mu_1)^p.
\end{align*}
Thus
\begin{align*}
W_p(\mu_s,\mu_t)\le (t-s)W_p(\mu_0,\mu_1).
\end{align*}
[guided]
We want to estimate the Wasserstein distance between the two intermediate measures $\mu_s$ and $\mu_t$. The definition of $W_p$ asks us to minimize over all couplings of $\mu_s$ and $\mu_t$, so it is enough to exhibit one admissible coupling and compute its cost.
The natural coupling is obtained by taking a random geodesic $\gamma$ distributed according to $\Pi$ and looking at its two positions $\gamma(s)$ and $\gamma(t)$. Formally, define
\begin{align*}
\pi_{s,t}:=(e_s,e_t)_\#\Pi.
\end{align*}
This is a probability measure on $X\times X$. Its first marginal is $(e_s)_\#\Pi=\mu_s$, and its second marginal is $(e_t)_\#\Pi=\mu_t$, so $\pi_{s,t}\in\Pi(\mu_s,\mu_t)$.
Because $W_p(\mu_s,\mu_t)^p$ is the infimum of $\int_{X\times X}d(x,y)^p\,d\rho(x,y)$ over all $\rho\in\Pi(\mu_s,\mu_t)$, the particular admissible plan $\pi_{s,t}$ gives the upper bound
\begin{align*}
W_p(\mu_s,\mu_t)^p\le \int_{X\times X} d(x,y)^p\,d\pi_{s,t}(x,y).
\end{align*}
Now we rewrite this integral using the definition of pushforward. Since $\pi_{s,t}=(e_s,e_t)_\#\Pi$, we have
\begin{align*}
\int_{X\times X} d(x,y)^p\,d\pi_{s,t}(x,y)=\int_{\operatorname{Geo}(X)} d(\gamma(s),\gamma(t))^p\,d\Pi(\gamma).
\end{align*}
The key point is that $\Pi$ is supported on constant-speed geodesics. Therefore, for every $\gamma\in\operatorname{Geo}(X)$ and every $0\le s\le t\le1$,
\begin{align*}
d(\gamma(s),\gamma(t))=(t-s)d(\gamma(0),\gamma(1)).
\end{align*}
Substituting this identity into the previous estimate gives
\begin{align*}
W_p(\mu_s,\mu_t)^p\le (t-s)^p\int_{\operatorname{Geo}(X)} d(\gamma(0),\gamma(1))^p\,d\Pi(\gamma).
\end{align*}
Finally, the assumption $(e_0,e_1)_\#\Pi=\pi$ says that the distribution of the pair of endpoints $(\gamma(0),\gamma(1))$ is exactly $\pi$. Hence
\begin{align*}
\int_{\operatorname{Geo}(X)} d(\gamma(0),\gamma(1))^p\,d\Pi(\gamma)=\int_{X\times X} d(x,y)^p\,d\pi(x,y).
\end{align*}
The coupling $\pi$ is optimal for $W_p(\mu_0,\mu_1)$, so
\begin{align*}
\int_{X\times X} d(x,y)^p\,d\pi(x,y)=W_p(\mu_0,\mu_1)^p.
\end{align*}
Combining these identities yields
\begin{align*}
W_p(\mu_s,\mu_t)^p\le (t-s)^pW_p(\mu_0,\mu_1)^p.
\end{align*}
Taking $p$th roots is valid because both sides are non-negative, and we obtain
\begin{align*}
W_p(\mu_s,\mu_t)\le (t-s)W_p(\mu_0,\mu_1).
\end{align*}
[/guided]
[/step]
[step:Force equality by applying the triangle inequality along the time chain]
Fix $0\le s\le t\le1$, and set
\begin{align*}
L:=W_p(\mu_0,\mu_1).
\end{align*}
Applying the previous step to the pairs $(0,s)$, $(s,t)$, and $(t,1)$ gives
\begin{align*}
W_p(\mu_0,\mu_s)\le sL.
\end{align*}
\begin{align*}
W_p(\mu_s,\mu_t)\le (t-s)L.
\end{align*}
\begin{align*}
W_p(\mu_t,\mu_1)\le (1-t)L.
\end{align*}
By the triangle inequality for the metric $W_p$ on $\mathcal P_p(X)$,
\begin{align*}
L=W_p(\mu_0,\mu_1)\le W_p(\mu_0,\mu_s)+W_p(\mu_s,\mu_t)+W_p(\mu_t,\mu_1).
\end{align*}
Combining this with the three upper bounds gives
\begin{align*}
L\le W_p(\mu_0,\mu_s)+W_p(\mu_s,\mu_t)+W_p(\mu_t,\mu_1)\le sL+(t-s)L+(1-t)L=L.
\end{align*}
Therefore both inequalities are equalities. Since each summand is bounded above by the corresponding summand in $sL+(t-s)L+(1-t)L$, equality of the total sum forces
\begin{align*}
W_p(\mu_s,\mu_t)=(t-s)L.
\end{align*}
Thus, for all $0\le s\le t\le1$,
\begin{align*}
W_p(\mu_s,\mu_t)=(t-s)W_p(\mu_0,\mu_1).
\end{align*}
If $s>t$, the same identity follows by symmetry of $W_p$. Hence
\begin{align*}
W_p(\mu_s,\mu_t)=|t-s|W_p(\mu_0,\mu_1)
\end{align*}
for all $s,t\in[0,1]$.
[/step]
[step:Conclude that the curve is a constant-speed Wasserstein geodesic]
Taking $s=0$ and $t=1$ in the preceding identity gives the prescribed endpoints:
\begin{align*}
\mu_0=(e_0)_\#\Pi
\end{align*}
and
\begin{align*}
\mu_1=(e_1)_\#\Pi.
\end{align*}
For arbitrary $s,t\in[0,1]$, the identity
\begin{align*}
W_p(\mu_s,\mu_t)=|t-s|W_p(\mu_0,\mu_1)
\end{align*}
says exactly that $[0,1]\ni t\mapsto\mu_t\in\mathcal P_p(X)$ has constant speed $W_p(\mu_0,\mu_1)$ with respect to the metric $W_p$. Therefore $(\mu_t)_{t\in[0,1]}$ is a constant-speed geodesic in $(\mathcal P_p(X),W_p)$ from $\mu_0$ to $\mu_1$.
[/step]
