Why is the distribution of mean firing rates assumed to be Gaussian? I'm reading "Principles of Neural Information Theory" by James V Stone and in section 3.5 he says that the distribution of firing rates (of a single neuron) is generally assumed to be approximately Gaussian. He proceeds to give a mathematical argument for this.

If we measure the activity during a period of $T$ seconds over intervals $\Delta t$, there are $N=T/\Delta t$ possible positions for spikes to occur. The probability of $n$ spikes occurring  this period is
$$p(n)=\frac{N!}{n!(N-n)!}P^nQ^{N-n}$$
Where $P$ is the probability fo a spike occurring and $Q$ is the probability of a spike not occuring. If $N$ and $T$ increase while firing rate $r$ stays constant, $p(n)$ approaches the Poisson distribution. If $T$ is held constant then a simple change of variables can be used to obtain a distribution for firing rate. If $N$ is large and $P$ is small, apparently this distribution is approximated by a Gaussian distribution.

I don't really understand this derivation. It would be nice if someone could provide an in-depth version of the argument with all steps full fleshed-out. I'm particularly confused about how one goes from p(n) to a Poisson distribution.
The distribution is not Gaussian, it's Poisson, as noted in your question. However, outside of certain edge cases (e.g. when there are very few spikes, or very few trials), a Poisson Distribution with rate parameter $\lambda$ looks very like a Gaussian with both mean and variance equal to $\lambda$ (standard deviation $= \sqrt{\lambda}$). This is useful, because if you treat the distribution as a Gaussian, you can apply methods such as the t-test or linear regression, without being too far wrong.

Here's some R code that shows the correspondance between the Poisson and the Gaussian.

library(tidyverse)
x = seq(0, 30, 1) # Possible spike counts
poisson.density = dopois(x, 10)
gaussian.density = dnorm(x, 10, sqrt(10))
d = data.frame(x = x, Poisson=poisson.density, Gaussian=gaussian.density) %&gt;%
  gather(Distribution, Density, -x)
ggplot(d, aes(x, Density, color=Distribution, linetype=Distribution)) +
  geom_path() +
  theme(legend.position=c(1,1), legend.justification=c(1,1)) +
  labs(x='Spike count')