Presumably you've solved this, but, in case not, it's because the PDE is a Kolmogorov backward equation, so the first and second order coefficients are the mean and variance of the underlying stochastic process being modeled.
In detail, consider a stochastic differential equation (which has a solution given by an Ito diffusion process):
$$ dp_t = \mu(p_t,t) dt + \sigma(p_t,t) dW_t $$
then the following system holds (under some conditions):
$$ -\frac{\partial}{\partial t} u(p,t) = \mu(p,t)\frac{\partial}{\partial p}u(p,t) + \frac{1}{2}\sigma^2(p,t)\frac{\partial^2}{\partial p^2} u(p,t) $$
where $u$ is the density of $p$ at $t$.
Note that the drift (infinitesimal mean) $M=\mu(p,t)$ and diffusion coefficient (infinitesimal variance) $V=\sigma^2(p,t)$ are as in the paper (except for the negative sign, which I assume is ignorable since he mostly seems to care only about the case when $\partial_t u\approx 0$ anyway).
Indeed, they are equivalently written:
\begin{align}
\mu(p,t) &= \lim_{\delta t\rightarrow 0}\frac{1}{\delta t} \mathbb{E}\left[ p_{t + \delta t} - p_t \mid p_t=p \right] =: M \\
\sigma^2(p,t) &= \lim_{\delta t\rightarrow 0}\frac{1}{\delta t} \mathbb{E}\left[ (p_{t + \delta t} - p_t)^2 \mid p_t=p \right] =: V \\
\end{align}
as Kimura writes.
Note that a useful approximation of the transition density is given by:
$$ \mathbb{P}[p_{t+\delta t} \mid p_t] \approx
\mathcal{N}(p_{t+\delta t}\mid p_t + \mu(p_t,t)\,\delta t, \sigma^2(p_t,t) \,\delta t) \tag{TD} $$
Ok, so everything above is just basic stochastic processes theory. Thus, if we have a stochastic model for the population dynamics, we can derive values for $M$ and $V$ from it (by computing its moments), and those will carry over to the backward Kolmogorov equation, upon which Kimura's work rests.
Here is where my ignorance of population dynamics shows. However, since Kimura mentions Fisher and Wright, I looked up the Wright-Fisher model.
It seems like Kimura is using the diffusion process approximation of the Wright-Fisher model.
This appears to be a well-studied and storied model that I cannot fully describe here; instead, I found the work by Tataru et al, Statistical Inference in the Wright–Fisher Model Using Allele Frequency Data to be an excellent description of it, though I do not pretend to understand much of it.
What's important, however, is that the change in genes (transition density) can be described by a binomial distribution. This can be approximated by a normal distribution:
$$ \mathbb{P}[p_{t+\delta t} \mid p_t] \approx \mathcal{N}(p_{t+\delta t}\mid p_t + a(p_t) \delta t,\, p_t (1-p_t) \delta t ) $$
using the standard approximation to the binomial.
This then gives us a forward Kolomogorov equation (not backward) written:
$$
\frac{\partial}{\partial t} u
=
- \frac{\partial}{\partial p}\left[a(p_t) u(p_t) \right]
+
\frac{1}{2}\frac{\partial^2}{\partial p^2} \left[ p_t(1-p_t)u(p_t) \right]
$$
This basically implies that $V=p(1-p)$.
( I noticed that another way to prove this is to notice that the Wright-Fisher approximate diffusion (without any selections, etc... so $a\equiv 0$) has an infinitesimal generator given by: $ \mathfrak{G} f(p) = p(1-p)\partial_{tt} f(p) / 2 $. This immediately implies $V=p(1-p)$. But may be less straightforward to understand. )
Confusingly, however, the paper has changed timescales (variables), so that $\delta t \leftarrow \Delta t / (2N)$, and then set $\delta t$ to $1$ (probably so that they would not have to write $2N$ everywhere).
If we undo this transform, we get
$$ \mathbb{P}[p_{t+\delta t} \mid p_t] \approx \mathcal{N}(p_{t+\delta t}\mid p_t + a(p_t) \delta t,\, p_t (1-p_t) \delta t/(2N) ) $$
If you compare this to our approximate transition density above (equation (TD)), you will see that this implies: $$ \sigma^2 = V = p(1-p)/[2N] $$
as desired.
Now, what is the infinitesimal mean, i.e. $a$ or $M$? Clearly this depends on the selection model, since it controls how the "environment" deterministically affects the process. Kimura describes this as a "constant selection advantage" with coefficient $s$.
The Tataru paper notes that the diffusion approximation to Wright-Fisher under genetic drift, mutation, and selection is given by:
$$
a(p) = -\nu p + \xi(1-p) + 2N\tau p(1-p)[h-(1-2h)p]
$$
If we (1) ignore mutation by setting $\nu=\xi=0$, (2) remove allelic dominance effects by setting $h=1/2$, and (3) define $s:= N\tau$, we get:
$$
a(p) = s p(1-p) =: M
$$
which of course we see by noting $M=a(p)$ matches $\mu$ in the equation (TD) above. (Note that the $2N$ transform also occurred here, but it was hidden inside $s$).
Thus, we have derived where Kimura's $M$ and $V$ come from, albeit probably not in the simplest possible manner.
All that remains is to derive the (steady-state) equation for $u$. I guess I'll do it for completeness.
Ignoring the steady-state subscripts, we get:
\begin{align}
G(x) &= \exp\left( -\int \frac{2M}{V} dx \right)
= \exp\left( -\int 4sN dx \right)
= \exp\left( -4sNx\right) \\[0.15cm]
u(p) &= \frac{\displaystyle \int_0^p G(x) dx}{\displaystyle \int_0^1 G(x) dx}
= \frac{\displaystyle \frac{1}{4Ns}\left[ \exp\left( -4sNx\right) \right]_0^p}{\displaystyle \frac{1}{4Ns}\left[ \exp\left( -4sNx\right)
\right]_0^1}
= \frac{\displaystyle -\left[ \exp\left( -4sNp\right) - 1 \right]}{\displaystyle -\left[ \exp\left( -4sN\right) -1
\right]} \\
&= \frac{1 - \exp(-4Nsp)}{1 - \exp(-4Ns)}
\end{align}
as required.
Apologies for any errors. (I am neither a population dynamics modeler nor a mathematician, so please point out any problems).
