Major update, based on your comment.
The code should be something like when using lme from nlme:
lme.model <- lme(dv ~ season * treat, random = ~1|rep, data=data)
where rep is a factor assigning unique codes to each of your 10 independent study sites (i.e., 5 per treatment), season is a factor indicating measurment quarter, and treatment is the treatment factor.
However, this will not give you a real ANOVA but a mixed model with one random effect (rep) and two fixed effects (season and between).
To fit a real ANOVA (namely one with one between- and one within-subjects factor, a so called split-plot design) you could use package afex:
require(afex)
anova <- ez.glm("rep", "dv", data = data, within = "season", between = "treat")
You could run this on each dv separately.
To analyze all dvs together you would need some multivariate analysis of which I am no expert.
Response prior to your comment:
If treat is your unit of observation (of which having two seems to be quite low), then the following code would be correct:
lme <- lme(dv ~ season, random = ~1|treat, data=data)
However, as said, having a repeated measures ANOVA with only two units of observation is pretty uncommon and seems a bad idea. If this is really your design (observed two treatments over several seasons), you are probably better off with other analyses, such as single-case analysis, e.g., here.
Best Answer
If I'm right in guessing that the other IV was not included in the ANOVA, then the most likely reason is that the two countries differed on the other IV, and that the countries only look 'significant' if the other IV is omitted. I wonder if the other IV is 'significant'? With sufficient confounding, it may not be. One last (but I think unlikely) possibility is that when you added the IV you lost 1 degree of freedom, and so if the IV were totally unrelated to the response, you would have lost a trivial amount of statistical power.