You can use glht for comparisons with a glm model object.

However, "tukey" in the mcp function simply sets up all pairwise contrasts for the grouping variable.  These are called Tukey contrasts.* **

I use glht for some things (http://rcompanion.org/rcompanion/h_01.html),

but I find lsmeans has so many advantages:

 - It can handle a variety of models and lists them in the documentation giving you some confidence it handles them appropriately.

 - It can handle interactions

 - It produces marginal-means (least square means)

 - It is easy to produce a compact letter display.

Attached is a txt file with an example, R code for using lsmeans with glm.

--------------

*  As an example of Tukey contrasts, see

 data(iris)

 library(multcomp)

 model = lm(Sepal.Length ~ Species, data=iris)

 glht(model, mcp(Species = "Tukey"))$linfct

**  At least for a one-way analysis of variance with a linear model, TukeyHSD, glht, and lsmeans can be made to produce the same results.

 if(!require(multcomp)){install.packages("multcomp")}

 if(!require(lsmeans)){install.packages("lsmeans")}

 if(!require(multcompView)){install.packages("multcompView")}

 data(iris)

 model = lm(Sepal.Length ~ Species, data=iris)

### Tukey HSD

 TukeyHSD(aov(model))

### glht

 library(multcomp)

 G = glht(model, mcp(Species = "Tukey"))

 confint(G)

### lsmeans

 library(lsmeans)

 leastsquare = lsmeans(model,

  pairwise ~ Species,

 adjust="tukey")

 confint(leastsquare$contrasts)

By the output seems like you try to perform a post hoc test on a continuous variable. For me this make no sense. Are you sure that your species_count_dry is categorical?. Another thing you should call TukeyHSD on a model that you created (i.e using aov or glm function) instead  on a data.frame object. Maybe the help tutorial in the packge multcomp could help you to understand better the TukeyHSD function.

The answer is even simpler.  It appears you can't use TukeyHSD with a glm object.

My next response has attached an example of using lsmeans with a glm object.

Note that estimates may not be in the original scale of the dependent variable. 

Read library(lsmeans); ?lsmeans .   Read ?lsmeans::models .

Note that the object GLM1 is a data.frame. So for me the  problem is not the use of TukeyHSD on the glm model. By the way TukeyHSD should works fine on his model since the way that him called glm is the same as perform a linear model and linear model is the same as aov model. 

I think the code below should work

require(multcomp)

glm_model

You can use glht for comparisons with a glm model object.

However, "tukey" in the mcp function simply sets up all pairwise contrasts for the grouping variable.  These are called Tukey contrasts.* **

I use glht for some things (http://rcompanion.org/rcompanion/h_01.html),

but I find lsmeans has so many advantages:

 - It can handle a variety of models and lists them in the documentation giving you some confidence it handles them appropriately.

 - It can handle interactions

 - It produces marginal-means (least square means)

 - It is easy to produce a compact letter display.

Attached is a txt file with an example, R code for using lsmeans with glm.

--------------

*  As an example of Tukey contrasts, see

 data(iris)

 library(multcomp)

 model = lm(Sepal.Length ~ Species, data=iris)

 glht(model, mcp(Species = "Tukey"))$linfct

**  At least for a one-way analysis of variance with a linear model, TukeyHSD, glht, and lsmeans can be made to produce the same results.

 if(!require(multcomp)){install.packages("multcomp")}

 if(!require(lsmeans)){install.packages("lsmeans")}

 if(!require(multcompView)){install.packages("multcompView")}

 data(iris)

 model = lm(Sepal.Length ~ Species, data=iris)

### Tukey HSD

 TukeyHSD(aov(model))

### glht

 library(multcomp)

 G = glht(model, mcp(Species = "Tukey"))

 confint(G)

### lsmeans

 library(lsmeans)

 leastsquare = lsmeans(model,

  pairwise ~ Species,

 adjust="tukey")

 confint(leastsquare$contrasts)

Thanks Salvatore S. Mangiafico for the guidance. It helped my analysis.

I needed a post hoc test for a continuous dependent variable and categorical independent variables. I had only two independent variables affecting the response variable additively. I could test the pairwise comparisons of levels with the two categorical variables. Moreover my response variable is non normal.

In the glm model object can we use still use the glht function for post hoc tests (Tukey contrasts) even if the dependent variable is non-normal?