One of the remarkable things about the Gaussian function is that its Fourier transform is also a Gaussian. A Gaussian blur means to convolve a function with the Gaussian. Convolution in the frequency (Fourier) domain is multiplication. Since the Gaussian goes to zero rapidly as you move away from zero, high-frequency components get attenuated whereas low-frequency stuff are preserved.
A pretty interesting point made here last year is that picking x, y, z coordinates from three orthogonal gaussian distributions results in a spherically symmetrical distribution. https://news.ycombinator.com/item?id=9446126
Yes, this is the separability property of the Gaussian. It makes it easy to compute Gaussian blur as compared to, say, the circular pillbox filter (lens blur or bokeh simulation) because you can just convolve it different times separately in each dimension. So, blurring an n by n image with a Gaussian of size h would only take O(n^2 h) instead of O(n^2 h^2). When h is large compared to log n, you can use the property I mentioned in my previous comment and do the Fast Fourier transform on the columns and rows of the image to improve time complexity to O(n^2 log n). It doesn't have to be a symmetrical distribution either; but in the general case with a multivariate Gaussian with some d-dimensional covariance matrix S, you'd have to rotate the data to align with the eigenvectors of S, which is lossy.
The rectangular pillbox filter also has the separable property. So does the parallelogram pillbox filter (albeit also rotated), which, by extension, allows us to do blurring by a hexagon (sum of three parallelograms) which can be used to simulate cool-looking bokeh [1].
[1] McIntosh, L., Bernhard E. Riecke, and Steve DiPaola. "Efficiently Simulating the Bokeh of Polygonal Apertures in a Post‐Process Depth of Field Shader." Computer Graphics Forum. Vol. 31. No. 6. Blackwell Publishing Ltd, 2012. http://ivizlab.sfu.ca/media/DiPaolaMcIntoshRiecke2012.pdf