I don't know if a mixed model is very appropriate (using the standard packages where the random effect structure is a linear predictor), unless you think the data at all time points should be exchangeable with each other in some sense (in which case the irregular intervals are a non-issue) - it wouldn't really be modeling the temporal autocorrelation in a reasonable way. It's possible you could trick lmer() into doing some sort of autogressive thing but how exactly you'd do that escapes me right now (I may not be thinking straight). Also, I'm not sure what the "grouping variable" would be that induces autocorrelation in the mixed model scenario.
If the temporal autocorrelation is a nuisance parameter and you don't expect it to be too large, then you could bin the data into epochs that are essentially disjoint from each other in terms of correlation (e.g. separate the time series at points where there are months of no data) and view those as independent replicates. You could then do something like an GEE on this modified data set where the "cluster" is defined by which epoch you are in, and the entries of the working correlation matrix are a function of how far apart the observations were made. If your regression function is correct, then you will still get consistent estimates of the regression coefficients, even if the correlation structure is misspecified. This would also allow you to model it as count data using, for example, the log-link (as one usually would in poisson regression). You could also build in some differential correlation between species, where each time point is viewed as a multivariate vector of species counts with some temporally decaying association between time points. This would require some pre-processing to trick the standard GEE packages into doing this.
If the temporal autocorrelation is not a nuisance parameter, I would try something more like a structured covariance model where you view the entire dataset as one observation of a big multivariate vector such that covariance between observations $Y_{s},Y_{t}$ on species $u,v$ is
$$ {\rm cov}(Y_{s}, Y_{t}) = f_{\theta}(s,t,u,v) $$
where $f$ is some parametric function known up to a finite number of parameters, $\theta$, along with a number of parameters to govern the mean structure. You might need to "build your own" for a model like this, but I'd also not be surprised if there are MPLUS packages to do things like this for count data.