having constructed a glm we can pass them through an anova as such:

anova(mnegbin1, mnegbin2, test = "Chisq")
Likelihood ratio tests of Negative Binomial Models

Response: count
                        Model    theta Resid. df    2 x log-lik.   Test    df      LR stat. Pr(Chi)
1 origin + substrate + sample 18.46502      1232       -2459.886                                   
2          substrate + sample 18.46502      1232       -2459.886 1 vs 2     0 -1.500666e-11       1

This gives us a comparison of the log likelihoods between the models, but can something similar be done for between residual deviances?:

summary(mnegbin1)
    Null deviance: 5686.16  on 1304  degrees of freedom
    Residual deviance:  890.24  on 1232  degrees of freedom

summary(mnegbin2)
    Null deviance: 5686.16  on 1304  degrees of freedom
    Residual deviance:  890.24  on 1232  degrees of freedom

In this instance they are of course identical but does that actually tell us anything?