[proofplan]
We prove the formula directly, without appealing to a general dynamical-plan theorem. First we check that the pushed-forward measures $\mu_t$ lie in $\mathcal{P}_p(\mathbb{R}^n)$ and have the prescribed endpoints. Then, for $0 \le s \le t \le 1$, we push the same optimal plan $\pi$ forward by $(T_s,T_t)$ to obtain a coupling of $\mu_s$ and $\mu_t$, whose transport cost is exactly $|t-s|^p$ times the original optimal cost. This gives the upper bound on $W_p(\mu_s,\mu_t)$, and the triangle inequality along the chain $\mu_0,\mu_s,\mu_t,\mu_1$ forces equality.
[/proofplan]
[step:Verify that the interpolated measures lie in $\mathcal{P}_p(\mathbb{R}^n)$ and have the correct endpoints]
For each $t \in [0,1]$, the map $T_t: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}^n$ is continuous, hence Borel measurable. Therefore the pushforward measure $\mu_t := (T_t)_\#\pi$ is a Borel probability measure on $\mathbb{R}^n$.
We now verify the finite $p$-moment condition. Since $t \in [0,1]$, the convexity inequality $|a+b|^p \le 2^{p-1}(|a|^p+|b|^p)$ gives, for every $(x,y) \in \mathbb{R}^n \times \mathbb{R}^n$,
\begin{align*}
|T_t(x,y)|^p \le 2^{p-1}\bigl((1-t)^p |x|^p + t^p |y|^p\bigr).
\end{align*}
Since $(\operatorname{pr}_1)_\#\pi=\mu_0$ and $(\operatorname{pr}_2)_\#\pi=\mu_1$, we obtain
\begin{align*}
\int_{\mathbb{R}^n} |z|^p \, d\mu_t(z) = \int_{\mathbb{R}^n \times \mathbb{R}^n} |T_t(x,y)|^p \, d\pi(x,y).
\end{align*}
Thus
\begin{align*}
\int_{\mathbb{R}^n} |z|^p \, d\mu_t(z) \le 2^{p-1}\int_{\mathbb{R}^n \times \mathbb{R}^n} \bigl(|x|^p+|y|^p\bigr) \, d\pi(x,y).
\end{align*}
Using the marginal identities again,
\begin{align*}
\int_{\mathbb{R}^n} |z|^p \, d\mu_t(z) \le 2^{p-1}\left(\int_{\mathbb{R}^n} |x|^p \, d\mu_0(x) + \int_{\mathbb{R}^n} |y|^p \, d\mu_1(y)\right) < \infty.
\end{align*}
Hence $\mu_t \in \mathcal{P}_p(\mathbb{R}^n)$.
At the endpoints, $T_0=\operatorname{pr}_1$ and $T_1=\operatorname{pr}_2$, so
\begin{align*}
\mu_0 = (T_0)_\#\pi
\end{align*}
and
\begin{align*}
\mu_1 = (T_1)_\#\pi.
\end{align*}
Thus the curve has the required endpoints.
[/step]
[step:Push the optimal plan forward to couple two intermediate measures]
Fix $s,t \in [0,1]$ with $s \le t$. Define the two-time interpolation map $S_{s,t}: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}^n \times \mathbb{R}^n$ by
\begin{align*}
S_{s,t}(x,y) = (T_s(x,y),T_t(x,y)).
\end{align*}
Let
\begin{align*}
\pi_{s,t} := (S_{s,t})_\#\pi.
\end{align*}
Then $\pi_{s,t}$ is a Borel probability measure on $\mathbb{R}^n \times \mathbb{R}^n$.
Let $\operatorname{pr}_1,\operatorname{pr}_2: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}^n$ denote the coordinate projections. Since $\operatorname{pr}_1 \circ S_{s,t}=T_s$ and $\operatorname{pr}_2 \circ S_{s,t}=T_t$, the functoriality of pushforward measures gives
\begin{align*}
(\operatorname{pr}_1)_\#\pi_{s,t} = (T_s)_\#\pi = \mu_s.
\end{align*}
Similarly,
\begin{align*}
(\operatorname{pr}_2)_\#\pi_{s,t} = (T_t)_\#\pi = \mu_t.
\end{align*}
Therefore $\pi_{s,t} \in \Pi(\mu_s,\mu_t)$.
[guided]
We need an admissible transport plan between $\mu_s$ and $\mu_t$. The construction should remember that both $\mu_s$ and $\mu_t$ were produced from the same pair $(x,y)$ sampled according to $\pi$. This is why we define a map that records both intermediate positions at once.
For fixed $s,t \in [0,1]$ with $s \le t$, define the map $S_{s,t}: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}^n \times \mathbb{R}^n$ by
\begin{align*}
S_{s,t}(x,y) = (T_s(x,y),T_t(x,y)).
\end{align*}
This map is Borel measurable because $T_s$ and $T_t$ are continuous maps from $\mathbb{R}^n \times \mathbb{R}^n$ to $\mathbb{R}^n$. We define
\begin{align*}
\pi_{s,t} := (S_{s,t})_\#\pi.
\end{align*}
Thus $\pi_{s,t}$ is the law of the pair of random points obtained by starting at $x$, ending at $y$, and observing the straight-line interpolation at times $s$ and $t$.
To prove that $\pi_{s,t}$ is a coupling of $\mu_s$ and $\mu_t$, we check its two marginals. Let $\operatorname{pr}_1,\operatorname{pr}_2: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}^n$ denote the coordinate projections. The first coordinate of $S_{s,t}(x,y)$ is $T_s(x,y)$, so
\begin{align*}
\operatorname{pr}_1 \circ S_{s,t} = T_s.
\end{align*}
Using the identity $(f \circ g)_\#\nu=f_\#(g_\#\nu)$ for pushforward measures, we get
\begin{align*}
(\operatorname{pr}_1)_\#\pi_{s,t} = (\operatorname{pr}_1)_\#(S_{s,t})_\#\pi = (\operatorname{pr}_1 \circ S_{s,t})_\#\pi = (T_s)_\#\pi = \mu_s.
\end{align*}
The same computation for the second coordinate gives
\begin{align*}
(\operatorname{pr}_2)_\#\pi_{s,t} = (\operatorname{pr}_2)_\#(S_{s,t})_\#\pi = (\operatorname{pr}_2 \circ S_{s,t})_\#\pi = (T_t)_\#\pi = \mu_t.
\end{align*}
Hence $\pi_{s,t}$ has first marginal $\mu_s$ and second marginal $\mu_t$, which is precisely the statement that $\pi_{s,t} \in \Pi(\mu_s,\mu_t)$.
[/guided]
[/step]
[step:Use the pushed-forward coupling to obtain the constant-speed upper bound]
Because $\pi_{s,t}$ is a coupling of $\mu_s$ and $\mu_t$, the definition of $W_p$ gives
\begin{align*}
W_p(\mu_s,\mu_t)^p \le \int_{\mathbb{R}^n \times \mathbb{R}^n} |u-v|^p \, d\pi_{s,t}(u,v).
\end{align*}
By the definition of $\pi_{s,t}$ as the pushforward of $\pi$ under $S_{s,t}$,
\begin{align*}
\int_{\mathbb{R}^n \times \mathbb{R}^n} |u-v|^p \, d\pi_{s,t}(u,v) = \int_{\mathbb{R}^n \times \mathbb{R}^n} |T_s(x,y)-T_t(x,y)|^p \, d\pi(x,y).
\end{align*}
For every $(x,y) \in \mathbb{R}^n \times \mathbb{R}^n$,
\begin{align*}
T_t(x,y)-T_s(x,y) = (t-s)(y-x).
\end{align*}
Therefore
\begin{align*}
|T_t(x,y)-T_s(x,y)|^p = |t-s|^p |x-y|^p.
\end{align*}
Combining this with the optimality of $\pi$,
\begin{align*}
W_p(\mu_s,\mu_t)^p \le |t-s|^p \int_{\mathbb{R}^n \times \mathbb{R}^n} |x-y|^p \, d\pi(x,y).
\end{align*}
Since $\pi$ is optimal,
\begin{align*}
W_p(\mu_s,\mu_t)^p \le |t-s|^p W_p(\mu_0,\mu_1)^p.
\end{align*}
Taking $p$-th roots gives
\begin{align*}
W_p(\mu_s,\mu_t) \le |t-s| W_p(\mu_0,\mu_1).
\end{align*}
[/step]
[step:Apply the triangle inequality to force equality]
For $0 \le s \le t \le 1$, apply the preceding upper bound to the pairs $(0,s)$, $(s,t)$, and $(t,1)$. This gives
\begin{align*}
W_p(\mu_0,\mu_s) \le s W_p(\mu_0,\mu_1).
\end{align*}
Also,
\begin{align*}
W_p(\mu_s,\mu_t) \le (t-s)W_p(\mu_0,\mu_1).
\end{align*}
And
\begin{align*}
W_p(\mu_t,\mu_1) \le (1-t)W_p(\mu_0,\mu_1).
\end{align*}
Since $W_p$ is the Wasserstein metric on $\mathcal{P}_p(\mathbb{R}^n)$, its triangle inequality gives
\begin{align*}
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*}
Substituting the three upper bounds yields
\begin{align*}
W_p(\mu_0,\mu_1) \le \bigl(s+(t-s)+(1-t)\bigr)W_p(\mu_0,\mu_1).
\end{align*}
Since $s+(t-s)+(1-t)=1$, every inequality in this chain must be an equality. In particular,
\begin{align*}
W_p(\mu_s,\mu_t) = (t-s)W_p(\mu_0,\mu_1).
\end{align*}
For arbitrary $s,t \in [0,1]$, exchanging $s$ and $t$ if necessary gives
\begin{align*}
W_p(\mu_s,\mu_t) = |t-s|W_p(\mu_0,\mu_1).
\end{align*}
Thus $(\mu_t)_{t \in [0,1]}$ is a constant-speed geodesic from $\mu_0$ to $\mu_1$, with interpolation formula $\mu_t=(T_t)_\#\pi=((1-t)\operatorname{pr}_1+t\operatorname{pr}_2)_\#\pi$.
[/step]