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).
One can also take an adelic viewpoint. When studying non-archimedean fields such as the p-adics $Q_p$, compact subgroups such as $Z_p$ play a pivotal role. On the reals, it seems the natural analogue of these compact subgroups are the Gaussians (cf. Tate's thesis). One can sort of justify the existence and central role of Gaussians on the grounds that the real number system "needs" something like the compact subgroups that its non-archimedean siblings enjoy, though this doesn't fully explain why Gaussians would then be exponentiated quadratic in nature.