Why they use KL divergence in Natural gradient? Natural gradient aims to do a steepest descent on the "function" space, a manifold that is independent from how the function is parameterized. It argues that the steepest descent on this function space is not the same as steepest descent on the parameter space. We should favor the former. 

Since, for example in a regression task, a neural net could be interpreted as a probability function (Gaussian with the output as mean and some constant variance), it is "natural" to form a distance on the manifold under the KL-divergence (and a Fisher information matrix as its metric). 

Now, if I want to be creative, I could use the same argument to use "square distance" between the outputs of the neural nets (distance of the means) which I think is not the same as the KL. 

Am I wrong, or it is just another legit way? Perhaps, not as good?
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.