There is a sense in which specifying an unknown probability distribution with
known mean and variance by a gaussian with the corresponding mean
and variance is the correct choice, and that is when those are the only things you are willing to say about the random variable. The gaussian distribution is the maximum entropy distribution with a given mean and variance.
It is the best probability distribution to use when you literally want to assume nothing else about your variable.
Gaussian Case
Let's say the mean and the variance are literally the only things you know about your variable. Then, we might be interested in which probability distribution best describes the variable. Surely we would want to make as general a choice as possible. This can be formulated. If we want to maximize our ignorance of the probability distribution subject to some constraints, we want to maximize the differential entropy of our distribution subject to those constraints.
In your case, this means we want to maximize:
$$ S = - \int dx \, p(x) \log p(x) $$
subject to the constraints
$$ \int dx \, p(x) = 1 \quad \int dx \, x \, p(x) = \mu \quad \int dx\, (x-\mu)^2 p(x) = \sigma^2 $$
The first constraint ensures that we have a proper probability density, and the second and third constraints are your observations about the mean and the variance. We can solve this with Lagrange multipliers:
$$ -\int dx \, p(x) \log p(x) + \lambda_0 \left[ \int dx\, p(x) - 1 \right] + \lambda_1 \left[ \int dx \, x \, p(x) - \mu \right] + \lambda_2 \left[ \int dx \, (x- \mu)^2 p(x) - \sigma^2 \right] $$
Where, taking the variation of $p(x)$ we obtain
$$ \int dx\, \left[ -\log p(x) - 1 + \lambda_0 + \lambda_1 x + \lambda_2 (x- \mu)^2 \right] = 0 $$
which is only satisfied if the integrand vanishes everywhere, suggesting:
$$ p(x) = \exp \left[ 1 + \lambda_0 + \lambda_1 x + \lambda_2 (x - \mu)^2 \right] $$
The only remaining problem is to determine the $\lambda$s from the constraints, from which we find
$$ p(x) = \frac{1}{\sqrt{ 2 \pi \sigma^2 } } \exp \left( - \frac{(x-\mu)^2 }{ 2\sigma^2 } \right) $$
The standard gaussian.
Generalization
What is really neat is that this generalizes. If you happen to have a set of functions $f_i$ for which you know the expectations under your probability distribution:
$$ \langle f_i (x) \rangle = \alpha_i $$
The probability distribution with the largest entropy consistent with those observations is:
$$ p(x) \propto \exp \left[\sum_i \lambda_i f_i(x) \right] $$
Where the proportionality constant and all of the $\lambda$s must be fixed to make all of the observed observations true. If you wanted a better distribution for your variable, you could go on to measure other quantities and reform your estimate for the probability distribution in this way.
Boltzmann Distribution
What is really interesting is that the probability distribution with maximum entropy that has as its only constraint a known expectation is the exponential distribution:
$$ p(x) \propto \exp \left( \frac{x}{\mu} \right) $$
which you might recognize as the Boltzmann distribution
This is not an accident. ET Jaynes would use this fact to build his principle of maximum entropy, and more generally formulate statistical mechanics in terms of information theory. For an nice introduction, consider his paper Information Theory and Statistical Mechanics [doi] [pdf]
