In general, under regularity conditions such as compactness, one can use standard minimax duality argument ($\min_y\max_x F(x,y)=\max_x\min_y F(x,y)$ -- for bilinear or, more generally, concave-convex functions $F$) to quickly show that
\begin{multline}
W(\mu_1, \mu_2, \dots, \mu_m):=\inf_\gamma\int c d\gamma \\
=\inf_\gamma^*\sup_{f_1,\dots,f_m}^*\Big(
\int c(x) \gamma(dx)
+\sum_{j=1}^m \int f_i d\mu_i-\sum_{j=1}^m \int f_i(x_i) \gamma(dx)
\Big) \\
=\sup_{f_1,\dots,f_m}^*\inf_\gamma^*\Big(
\int c(x) \gamma(dx)
+\sum_{j=1}^m \int f_i d\mu_i-\sum_{j=1}^m \int f_i(x_i) \gamma(dx)
\Big) \\
=\sup
\sum_{j=1}^m \int f_i d\mu_i,
\end{multline}
where $\inf_\gamma$ is taken over all probability measures $\gamma$ with marginals $\mu_1, \mu_2, \dots, \mu_m$, $\inf\limits_\gamma^*$ is taken over all (nonnegative) measures $\gamma$ on $\mathcal{X}^m$,
$\sup\limits_{f_1,\dots,f_m}^*$ is taken over all appropriately integrable real-valued measurable functions $f_1,\dots,f_m$ on $\mathcal{X}$,
$\sup$ is taken over all such functions $f_1,\dots,f_m$ on $\mathcal{X}$ that satisfy the condition
\begin{equation}
c(x)\ge \sum_{j=1}^m f_i(x_i)\text{ for all }x\in\mathcal{X}^m,
\end{equation}
$x=(x_1,\dots,x_m)$, $dx=dx_1\times\dots\times dx_m$.
The second equality above follows because
\begin{multline}
\sup_{f_1,\dots,f_m}^*\Big(
\sum_{j=1}^m \int f_i d\mu_i-\sum_{j=1}^m \int f_i(x_i) \gamma(dx)
\Big) \\
=\begin{cases}
0&\text{ if the marginals of $\gamma$ are $\mu_1, \mu_2, \dots, \mu_m$}, \\
\infty&\text{ otherwise.}
\end{cases}
\end{multline}
The third equality above is the minimax duality.
The last equality there
follows because
\begin{multline}
\inf_\gamma^*\Big(
\int c(x) \gamma(dx)
-\sum_{j=1}^m \int f_i(x_i) \gamma(dx)
\Big) \\
=\begin{cases}
0&\text{ if $c(x)\ge \sum_{j=1}^m f_i(x_i)\text{ for all }x\in\mathcal{X}^m$}, \\
-\infty&\text{ otherwise.}
\end{cases}
\end{multline}
So, in general $\sum_{j=1}^m \int f_i d\mu_i$ and the condition $c(x)\ge \sum_{j=1}^m f_i(x_i)$ for all $x\in\mathcal{X}^m$ replace, respectively,
$\mathbb{E}_{\mu}f(X)-\mathbb{E}_\nu f(X)=\int f d\mu-\int f d\nu$
and the condition $|f(x)-f(y)|\leq d(x, y)$ for all $x, y\in \mathcal{X}$
in the definition of the Wasserstein distance.
Even in the case $m=2$, the reduction from two functions $f_1$ and $f_2$ in the condition $c(x)\ge \sum_{j=1}^2 f_i(x_i)$ for all $x\in\mathcal{X}^2$ to one function $f$ in the condition $|f(x)-f(y)|\leq d(x, y)$ for all $x, y\in \mathcal{X}$ is based on the special properties of a metric, available only if $c=d$, a metric. This is how it is done: Suppose $d(x_1,x_2)\ge \sum_{j=1}^2 f_i(x_i)$ for all $x\in\mathcal{X}^2$, that is, $f(s):=\inf\limits_t[d(s,t)-f_2(t)]\ge f_1(s)$ for all $s$. Also, $f(s)\le d(s,s)-f_2(s)=-f_2(s)$ for all $s$, so that $f_1\le f$ and $f_2\le-f$ and hence
\begin{equation}
\sum_{j=1}^2 \int f_i\, d\mu_i\le \int f\, d\mu_1-\int f\, d\mu_2.
\end{equation}
Also,
$$d(s,t)+f(t)=\inf\limits_u[d(s,t)+d(t,u)-f_2(u)]
\ge\inf\limits_u[d(s,u)-f_2(u)]=f(s)$$
and hence $d(s,t)\ge f(s)-f(t)$ for all $s,t$ or, equivalently,
$|f(t)-f(s)|\le d(s,t)$ for all $s,t$.
