At first, you must reorganize your table:
Months Replicate2 Species.1 Species.2 Species.3 Species.4
Jan-12 BranA 778 41 19 1
Jan-12 BranB 800 18 18 2
Jan-12 BranC 537 48 20 2
Feb-12 BranA 465 29 19 1
Feb-12 BranB 444 45 14 1
Feb-12 BranC 671 13 10 2
...
# reorganisation procedure
# (need to be optimalized according to your full datatable!!!)
# create vectors...
my.data <- read.table("spec.csv", header = TRUE, sep = ",") # load your data
spec_abund <- c(my.data$Species.1, my.data$Species.2,
my.data$Species.3, my.data$Species.4) # species abundance
sites <- rep(my.data$Replicate2, times = 4) # sites
my.date <- rep(my.data$Months, times = 4) # date of collection
spec.no <- rep(c("spec1","spec2","spec3","spec4"),
each = nrow(my.data)) # spec.ID
# ...and bind them into new data.frame
my.reorg <- data.frame(abundance = spec_abund, sites = sites,
months = my.date, spec.no = spec.no) # final table
It should looks like this:
my.reorg
abundance sites months spec.no
1 778 BranA Jan-12 spec1
2 800 BranB Jan-12 spec1
3 537 BranC Jan-12 spec1
4 465 BranA Feb-12 spec1
5 444 BranB Feb-12 spec1
6 671 BranC Feb-12 spec1
...
Than try to fit the model:
aov.out <- aov(abundance ~ sites * months + Error(spec.no/months), data=my.reorg)
summary(aov.out)
It looks like that sites are non-significantly different (in species composition) from each other during the "whole time".
Is this what you need?