Quadratic (or bilinear) forms appear naturally throughout mathematics, for instance via inner product structures, or via dualisation of a linear transformation, or via Taylor expansion around the linearisation of a nonlinear operator.  The Laplace-Beltrami operator and similar second-order operators can be viewed as differential quadratic forms, for instance.

A Gaussian is basically the multiplicative or exponentiated version of a quadratic form, so it is quite natural that it comes up in multiplicative contexts, especially on spaces (such as Euclidean space) in which a natural bilinear or quadratic structure is already present.

Perhaps the one minor miracle, though, is that the Fourier transform of a Gaussian is again a Gaussian, although once one realises that the Fourier kernel is also an exponentiated bilinear form, this is not so surprising.  But it does amplify the previous paragraph: thanks to Fourier duality, Gaussians not only come up in the context of spatial multiplication, but also frequency multiplication (e.g. convolutions, and hence CLT, or heat kernels).