I have done an experiment of how different animal species affect nutrient fluxes in sediments. I had a number of experimental units (sediment boxes) to which different animal species were added. I then measured nutrient fluxes in these units. Each unit was measured five times. I have used the ezANOVA package for R to do a repeated measures ANOVA with two factors (between subjects factor “Species” and within subject factor “Time”). 
Here's an example of my input:

 > ezANOVA(data=NoP_3_5,dv=.(AcPO4),wid=.(Subject),within=.(Time),between=.(Species),return_aov=T)
 $ANOVA
 Effect DFn DFd F p p<.05 ges
 2 Species 2 6 10.60384830 1.072453e-02 * 0.50875722
 3 Time 2 12 27.88590570 3.081718e-05 * 0.76667542
 4 Species:Time 4 12 0.08093179 9.867291e-01 0.01871588

 $`Mauchly's Test for Sphericity`
 Effect W p p<.05
 3 Time 0.1439855 0.007866779 *
 4 Species:Time 0.1439855 0.007866779 *

 $`Sphericity Corrections`
 Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
 3 Time 0.5387889 0.001349263 * 0.5630403 0.001103048 *
 4 Species:Time 0.5387889 0.933492515 0.5630403 0.939138882 

 $aov

 Call:
 aov(formula = formula(aov_formula), data = data)

 Grand Mean: -0.004276663 

 Stratum 1: Subject

 Terms:
 Species Residuals
 Sum of Squares 0.0004410742 0.0001247870
 Deg. of Freedom 2 6

 Residual standard error: 0.004560464 
 Estimated effects may be unbalanced

 Stratum 2: Subject:Time

 Terms:
 Time Species:Time Residuals
 Sum of Squares 0.0013994206 0.0000081229 0.0003011028
 Deg. of Freedom 2 4 12

 Residual standard error: 0.005009181 
 Estimated effects may be unbalanced
 

I want to use the $aov output to do a Tukey HSD post hoc test of the between factor ("Species"). However, it does not work:


 mod <- ezANOVA(data=NoP_3_5,dv=.(AcPO4),wid=.(Subject),within=.(Time),between=.(Species),return_aov=T)
 > TukeyHSD(mod$aov)
 Error in UseMethod("TukeyHSD") : 
 no applicable method for 'TukeyHSD' applied to an object of class "list"

I have also tried to use the aov-command directly:

 > mod2 <- aov(AcPO4~(Species)+Error(Subject/Time)+(Species),data=NoP_3_5)
 > TukeyHSD(mod2)
 Error in UseMethod("TukeyHSD") : 
 no applicable method for 'TukeyHSD' applied to an object of class "c('aovlist', 'listof')"


I have tried to use lme() to specify the anova-model with the aim to produce something that TukeyHSD() without success.

Any help would be greatly appreciated!