Let's consider what these statistics are, and then get into what mean substitution does.
Definitions
Ignoring the missing values for a moment, you have a dataset of numbers $x_1, x_2, \ldots, x_n,$ not all of them equal. The SD $\sigma,$ skewness $\beta,$ and kurtosis $\kappa$ of this dataset all depend on the residuals relative to the mean
$$m = (x_1 + x_2 + \cdots + x_n)/n.$$
These residuals are
$$y_i = x_i - m.$$
They sum to zero (because the $x_i$ and the $n$ subtracted copies of $m$ exactly balance). In these terms there are various formulas for things called "standard deviation" (SD), "skewness," and "kurtosis." Because the details don't matter, let's express these formulas as
$$\mu_k = \frac{1}{C_{n,k}} \left(y_1^k + y_2^k + \cdots + y_n^k\right)$$
where $k=2$ gives the square of the SD (that is, the variance), $k=3$ gives the third central moment, and $k=4$ gives the fourth central moment. The denominators $C_{n,k}$ increase with $n,$ are close to $n,$ and (on a relative basis) approach $n$ as $n$ grows large. See https://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm for a rundown of some of the formulas to see the details I'm avoiding here.
Finally, the statistics in question are
$$SD = \sigma= \sqrt{\mu_2};\quad \beta = \frac{\mu_3}{\sigma^3} = \frac{\mu_3}{\mu_2^{3/2}};\quad \kappa = \frac{\mu_4}{\sigma^4}= \frac{\mu_4}{\mu_2^{2}}.$$
Analysis
Consider those missing values and suppose there are $a$ of them. "Mean substitution" introduces additional "data" $x_{n+1}=x_{n+2}=\cdots=x_{n+1}=m$ all equal to the original mean. Obviously the new mean is the same as the old, which implies all the residuals are the same. The new residuals are $y_{n+i} = x_{n+i}-m = m-m= 0.$ Consequently, mean substitution does not change any of the sums of powers in the definitions of the central moments $\mu_k.$
However, those pesky denominators do change: where originally $C_{n,k} \approx n$ is used, with the mean-substituted data $C_{n+a,k}\approx n+a$ will be used in its place. As a result,
Mean-substitution of $a$ missing values multiplies the original central moment $\mu_k$ by $C_{n,k}/C_{n+a,k}.$
Notice this implies a decrease in the size of every central moment. This is the crux of the matter.
Results
Plugging the preceding result into the formulas for the statistics and using primes for the new values, we see
$$\begin{aligned}
\sigma^\prime &= \sqrt{\frac{C_{n,2}}{C_{n+a,2}}}\,\sigma &\approx \sqrt{\frac{n}{n+a}}\,\sigma &\lt \sigma;\\
\beta^\prime &= \frac{C_{n,3}}{C_{n+a,3}}\left(\frac{C_{n+a,2}}{C_{n,2}}\right)^{3/2} \beta &\approx \sqrt{\frac{n+a}{n}}\beta;\\
\kappa^\prime &= \frac{C_{n,4}}{C_{n+a,4}}\left(\frac{C_{n+a,2}}{C_{n,2}}\right)^2 \kappa &\approx \frac{n+a}{n}\kappa & \gt \kappa.
\end{aligned}$$
That is,
The SD $\sigma$ decreases because a bunch of zero residuals have been thrown into the mix;
So the skewness (if not originally zero) increases because the decrease in $\sigma^3$ in the denominator overwhelms the decrease in $\mu_3$ from the inclusions of (cubes) of those zero residuals;
And the kurtosis increases even more from the decrease in $\sigma^4.$
By "increase" in skewness I mean an increase in its size: its sign remains the same.
Notice that these are universal patterns on a relative scale: that is, we can predict exactly how $\sigma^\prime/\sigma,$ etc. will behave as a function of $n$ and $a$ once we have selected our denominators $C_{n,k}$ for computing these statistics. With $C_{n,k}=n$ in all cases (these are the empirical formulas) and $n=4$ non-missing values, here is what these plots look like (with positive skewness):
For those who like to experiment, here is the R code to produce the illustration. It is written to parallel the mathematical exposition here rather than for efficiency.
mu <- function(x, k) mean((x - mean(x))^k) # Central moment of `x` of order `k`
moment.std <- function(x, k) mu(x, k) / mu(x, 2)^(k/2) # Standardized central moment of `x`
sigma <- function(x) sqrt(mu(x, 2))
beta <- function(x) moment.std(x, 3)
kappa <- function(x) moment.std(x, 4)
x <- c(0,1,2,5) # Any set of values of nonzero skewness will serve as an illustration
a <- 0:10 # Number of missing values to examine
stats <- sapply(a, function(a) {
z <- c(x, rep(mean(x), a)) # Mean substitution of `a` values
c(Sigma=sigma(z), Beta=beta(z), Kappa=kappa(z)) # The three statistics
})
stats <- stats/stats[, 1] # Relative values of the statistics
par(mfrow=c(1,3))
plot(a, stats["Sigma",], ylab="", main=expression(sigma * "'" /sigma))
plot(a, stats["Beta",], ylab="", main=expression(beta * "'" /beta))
plot(a, stats["Kappa",], ylab="", main=expression(kappa * "'" /kappa))
par(mfrow=c(1,1))
