How to do MANCOVA in R?

I have two groups of persons, GRP0 and GRP1, on which I measured three continuous variables: VAR1, VAR2 and VAR3.

I would like to use Mancova in R with: - VAR1, VAR2 and VAR3 as outcome variables - GRP={0,1} as predictor variable - age and gender as covariates

What would be the correct way to formulate this model in R?

Also the measured were carried out on 100 instances (which are serially correlated) so any help on how to apply permutation-based multiple-comparison correction on top would be ideal.

Thanks a lot

Y <- cbind(VAR1, VAR2, VAR3)
fit <- manova(Y ~ GRP+age+gender)
summary(fit, test="Pillai")


Other test options are "Wilks", "Hotelling-Lawley", and "Roy".

summary.aov(fit) # for univariate statistics.