The problem of Gaussian approximation for the wireless multi-access interference distribution in large spatial wireless networks is addressed. First, a principled methodology is presented to establish rates of convergence of the multi-access interference distribution to a Gaussian distribution for general bounded and power-law decaying path-loss functions. The model is general enough to also include various random wireless channel dynamics such as fading and shadowing arising from multi-path propagation and obstacles existing in the communication environment. It is shown that the wireless multi-access interference distribution converges to the Gaussian distribution with the same mean and variance at a rate 1√λ, where λ > 0 is a parameter controlling the intensity of the planar (possibly non-stationary) Poisson point process generating node locations. An explicit expression for the scaling coefficient is obtained as a function of fading statistics and the path-loss function. Second, an extensive numerical and simulation study is performed to illustrate the accuracy of the derived Gaussian approximation bounds. A good statistical fit between the interference distribution and its Gaussian approximation is observed for moderate to high values of λ.