I am working on an RNA-seq dataset for differential expression: 4 environments, 4 biological replicates per environment. There are also 4 batches, but they are orthogonal to environments, so all environments are equally represented among all batches.

I lean towards EdgeR. I find it more user friendly and have used it more often, easier to select for a common-sense effect size in addition to p value, it has a robust way of estimating variation, uses negative binomial distribution (but this is experimental in lme4) and so on.

However, some senior colleagues wants to use the mixed model approach of lme4 so that batch will be a random factor. Why does EdgeR not allow the use of random factors? What kind of arguments could I use in favor of EdgeR over lme4? Or if that shows my bias too much, what are the strengths of each approach?

From what I understand, mixed effect models are HARD; hard to fit and hard to get inferences from. I once looked at lme4's internal code and I came away with a headache. It certainly wouldn't be trivial to add support for random effects to edgeR. In any case, here are a couple of points you can present to your colleagues:

• You don't need random effects. The difference between fixed and random effects is less to do about the expectation of whether something is random or not (I mean, what isn't random when you do an experiment?) but rather about whether or not you need to share information across levels of the random effect (in this case, batches). This avoids the need to fully estimate a coefficient for each batch, as the estimates are shrunk towards the distribution from which they were assumed to be sampled. This preserves information for use elsewhere, which is important, e.g., when your batches are confounded with your variable of interest. Alternatively, if you want to know stuff about the underlying distribution, then random effects can also be useful. But neither of these things are an issue for your data set - the batches are orthogonal to the environment, and you don't really care about the batches other than to get them out of the way. From edgeR's side, type I error control will still be maintained for the environment comparisons when you model the batches as fixed effects.
• The use of an underlying common distribution for the random effect can run into problems when one batch is an outlier or if you have some further systematic differences between your batches (e.g., some of the batches were processed differently from the others). This will interfere with estimation of the variance of the random effect, whereas the fixed effect won't care.
• edgeR shares information across genes to improve the stability of the dispersion estimates, which increases the power to detect DE. This is not performed in lme4 - obviously, because that package wasn't designed for genome-wide studies. Any gain in power you might get from sharing information to estimate the batch coefficients in lme4 would probably be balanced out (and then some) by the gain in power in edgeR from increased precision of the dispersion estimates. In fact, I'm not sure that the variance of the random effect is estimated reliably when you have few batches; my personal tests with 6 batches resulted in loss of type I error control with glmer, and I suspect if you only have 4 batches, you'd probably get something similar.
• edgeR is faster than running a loop over all genes to fit a GLMM to each. I've tried it, and it took hours.

If you must use random effects for RNA-seq data, have a look at duplicateCorrelation with voom and limma. This is fast, shares information across genes to improve estimation, and generally does the job pretty well.

Well, there are a lots of reasons why using lme4 for RNA-seq data would be unstable and unreliable (see Aaron's items 2 - 4).

For your experiment, though, there is an additional reason to avoid lme4. There is actually no theoretical advantage for using random effects for a fully balanced design such as yours (with batches orthogonal to environments). You will get essentially full information about the environmental effects from the edgeR model with batch as a fixed effect. The theoretical reason for using a random effect model is to recover inter-batch information about the treatments but, for the fully balanced design, there is no inter-batch information to recover. (Aaron has said essentially the same thing in his item 1.)

If you were analyzing a univariate normal response variable, then the fixed and random effects models would give identical results for your design. For glm mixed models, everything is approximate, but the same general principles apply.