If we consider just one mode and we compute the contribution from this mode to the $n$-th moment $\langle\mu_n\rangle = \langle E(0)^n\rangle$ (the state is the vacuum) and we also compute the moments of a Gaussian distribution (zero mean and unit variance) we observe that the values are the same.
Do you think that this is a conclusive proof for saying that the vacuum fluctuations are Gaussian?
As innisfree points out in a comment, a distribution is gaussian if and only if all its moments are gaussian. In the case of a free bosonic field, Wick's theorem guarantees that all the $n$-point functions are gaussian, as can be seen by comparing the theorem to Isserlis' theorem. In the end it all boils down to the fact that $n$-point functions of free fields are formally gaussian moments[1]: $$ G_n\propto \int\mathrm d\varphi\ \varphi_1\varphi_2\cdots\varphi_n\ \mathrm e^{iS[\varphi]} $$ where $S$ is quadratic, and $\mathrm d\varphi$ is the measure over the space of field configurations (see Functional integration).
But there is a very important point to be made: quantum field theory is not statistical mechanics, regardless of the obvious formal analogies. The correlation functions of QFT do not measure correlations in the statistical sense of the word; the $G_n$ above are not moments of a distribution.
Moreover, as Mark Mitchison mentions in another comment above, the real word is not free and thus the true correlation functions are not gaussian. The action $S$ is not quadratic in general. And the "true configuration" of the universe is not the vacuum state. Therefore, the discussion above does not really apply to the real world.
[1] to make the analogy more clear one should Wick rotate the path integral.
