(1) I adjusted the claim to "..$x \notin \Gamma$ ...". The whole argument assumes $x \ne x_\Gamma$. (2) Note $f$ is $1$-Lipschitz and for general $1$-Lipschitz functions it holds that the gradient is equal to $e$ whenever a partial in direction $e$ is equal to $1$. Alternatively, observe that on the one hand $f(y)- f(x)$= -\|y-x\|$ by choice of $y$. On the other hand $f(y) = f(x) + v\cdot (x-y) + o(\|x-y\|) by differentiability in $x$. So if $y$ is close to $x$ then necessarily $v = e$ by the "equality case" of Cauchy-Schwarz.

In the answer I added that phi is integer-valued. This follows from formula (5.17) in Villani’s book and the assumption that c is integer-valued. I don’t claim that a given dual solution is integer-valued. Most likely it is possible to glue together two solutions. One that is constant and another one that has exactly oscillation equal to 1. In a case i could imagine (I don’t have time to construct such an example), there might be a dual solution that has values 0,0.5 and 1.

I don’t understand the problem as the argument does not claim any sharpness. Maybe it helps to clarify the notation: the constant function equal to 0 has values in $\{0,1\}$. The example given in your comment just verifies that $\varphi$ has values in $\{-1,0\}$.

Ah, this follows from the construction that the sum and difference of c evaluated at different couples is an integer. Thus phi is integer-values.

Was the right hand side. Was just a typo. This is true as c has values 0 or 1 so the difference can have at most the mentioned three values. The inequality then shows that the oscillation of phi is at most 1. Knowing phi at x0 is 0 gives the claims.

I changed to a direct argument.

Since the value question was "already" answered in the question itself, I didn't include it. Now there is a short statement.