Intuitively, why is the Gaussian the Fourier transform of itself?

The generalization of this phenomenon, from a probabilistic standpoint, is the Wiener-Askey Polynomial Chaos.

In general, there is a connection between orthogonal polynomial families in the Askey scheme and probability distribution/mass functions.

Orthogonality of these polynomials can be shown in an inner product space using a weighting function -- a weight function that typically happens to be, within a scale factor, the pdf/pmf of some distribution.

In other words, we can use these orthogonal polynomials as a basis for a series expansion of a random variable:

$$z = \sum_{i=0}^\infty z_i \Phi_i(\zeta).$$

The random variable $\zeta$ belongs to a distribution we choose, and the orthogonal polynomial family to which $\Phi$ belongs follows from this choice.

The deterministic coefficients $z_i$ can be computed easily by using Galerkin's method.

So, yes. There is a very deep connection in this regard, and it is extremely powerful, particularly in engineering applications. Strangely, many mathematicians do not know this relationship!

See also: http://www.dtic.mil/cgi-bin/GetTRDoc?AD=ADA460654 and the Cameron-Martin Theorem.

There's a simple reason why taking a Fourier-like transform of a Gaussian-like function yields another Gaussian-like function. Consider the property $$\mathcal{T}[f^\prime](\xi) \propto \xi \hat{f}(\xi)$$ of a transform $\mathcal{T}$. We will call an invertible transform $\mathcal{F}$ "Fourier-like" if both it and its inverse have this property.

Define a "Gaussian-like" function as one with the form $$f(x) = A e^{a x^2}.$$ Functions with this form satisfy $$f^\prime(x) \propto x f(x).$$ Taking a Fourier-like transform of each each side yields $$\xi \hat{f}(\xi) \propto \hat{f}^\prime(\xi).$$ This is has the same form as the previous equation, so it is not surprising that its solutions have the Gaussian-like form $$\hat{f}(\xi) = B e^{b \xi^2}.$$