Why they use KL divergence in Natural gradient?
The KL divergence has slightly different interpretations depending on the context. The related Wikipedia article contains a section dedicated to these interpretations. Independently of the interpretation, the KL divergence is always defined as a specific function of the cross-entropy (which you should be familiar with before attempting to understand the KL divergence) between two distributions (in this case, probability mass functions)

\begin{align}
D_\text{KL}(P\parallel Q) 
&amp;= -\sum_{x\in\mathcal{X}} p(x) \log q(x) + \sum_{x\in\mathcal{X}} p(x) \log p(x) \\
 &amp;= H(P, Q) - H(P)
\end{align}
where $H(P, Q)$ is the cross-entropy of the distribution $P$ and $Q$ and $H(P) = H(P, P)$.

The KL is not a metric, given that it does not obey the triangle inequality. In other words, in general, $D_\text{KL}(P\parallel Q) \neq D_\text{KL}(Q\parallel P) $.

Given that a neural network is trained to output the mean (which can be a scalar or a vector) and the variance (which can be a scalar, a vector or a matrix), why don't we use a metric like the MSE to compare means and variances? When you use the KL divergence, you don't want to compare just numbers (or matrices), but probability distributions (more precisely, probability densities or mass functions), so you will not compare just the mean and the variance of two different distributions, but you will actually compare the distributions. See the example of the application of the KL divergence in the related Wikipedia article.