Because the resulting model is useful. Sure, you could fit an elaborate model that includes all theoretical contributions to experimental and biological variability. Or, you could just fit a simple NB-based model that accounts for the major features of the data, i.e., discrete and overdispersed counts. See George Box's quote that "All models are wrong but some are useful". Will the expression levels of a gene be Gamma-distributed in reality? Almost definitely not, but for a continuous, unimodal (in bulk RNA-seq data) distribution of positive expression values, the Gamma distribution is close enough.

In short, don't think of it as a streetlight at night. Instead, it's more like a sunny day on a meadow with with a couple of trees. While some of the field is in shadow (corresponding to special cases that aren't handled well by the NB model, e.g., bimodal distribution of expression values in single-cell RNA-seq data), most of it can still be searched, even when drunk.

P.S. Alternative Poisson mixtures have been studied for other purposes, e.g., Poisson-Beta models for transcriptional bursting, Poisson-log-normal models for network modelling. But they involve more effort to fit and get inferences from, and are mainly used because of some property of the mixing distribution that is relevant to the specific application (e.g., on/off rates as parameters of the Beta distribution). For routine analyses, if something simpler does the job, why not use it?

As far as I know, the analytic/computational convenience is the primary justification. We're talking about modeling the biological variation here, and in theory the distribution of biological variation could be anything. It depends on the biological system in question. Unfortunately, most RNA-seq experiments have very few replicates, so analysis isn't possible without making some kind of parametric assumption about the biological variation. With many more replicates, you could use a non-parametric approach and avoid such assumptions, but with the typical number of replicates, this approach has very low power.

If you don't like the assumption of a gamma distribution, there are packages that offer more flexible distributions, such as the tweeDESeq, which adds a third parameter to the NB distribution that I believe allows fatter or thinner tails on the distribution biological variation than the gamma can accommodate. Again, though, estimating more parameters requires more replicates. There is also limma's voom function, which ignores the principle of fitting a mixture of Poisson distributions and simply fits a linear model (i.e. using the normal distribution) with observation weights derived from the empirical mean-variance trend. In both voom and in the case of edgeR and others choosing the gamma distribution to model the biological variation, the key idea is that the correctness of an inference is far more sensitive to the magnitude of the variance than the shape of the distribution. This is borne out by numerous studies in which the humble t-test is shown to be quite robust even in the face of substantial deviations from normality in the underlying data, and is also supported by the fact that limma-voom and edgeR's QL method give quite similar p-values on a good quality data set, despite the substantial differences in parametric assumptions.

Regarding the idea that limiting oneself to the gamma distribution represents an unreasonable limiting of once's search space, we can ask in what cases the gamma distribution will fail to capture the true magnitude of variability. I would assume that NB GLM methods would overestimate the significance (i.e. underestimate the FDR) when the biological variation has a heavy-tailed distribution and the gamma cannot fit those fat tails, and vice versa for light-tailed distributions. Additionally, there is probably room for error on distributions that are heavily skewed in one direction or the other.

@Aaron, @Ryan: I don't have an alternative to the GD that I think would be better.  My question was simply: out all the possible distributions one could have picked to model the uncertainty in the Poisson parameter, why the GD?  After reading more on the matter, I think the answer does boil down to analytic convenience, the GD being the conjugate prior of the Poisson, but in a Bayesian context this is the natural choice.  IOW, in the context of Bayesian updating, many weighting distribution would be comparably adequate; one may as well choose the one that is analytically most convenient, and that's the GD in this case.