The negative binomial has two parameter in its distribution. Neg bin has a scale and a probability parameter. I’d imagine the scale parameter estimated in poisson regression is only one of them. Does running poisson regression give me those two parameters in the negbin?
It is indeed the case that if you use family=negative.binomial (from the MASS package) that you need to specify theta. You could try running GEE fits in a loop or an optimizer with a range of theta values and selecting the one that maximizes some goodness-of-fit statistic (this is what MASS::glm.nb does).
I would suggest looking into Zuur et al's book on mixed models in ecology to see if there is anything useful there.