Nested test to remedy pseudoreplication

#1
Hello and thanks for taking the time to read this.

Data structure:
I have two datasets from two protected areas that differ in protection status. Both areas contain 43 and 37 sites each.

Question:
I would like to know which test would be the best for testing whether the PA status has had an effect on 1) the first axis of a PCoA (principal coordinates analysis) - i.e. species composition turnover (derived by constructing a bray curtis dissimilarity matrix) and 2) species richness per site (a continuous variable).

Problem:
I understand that there is pseudoreplication present in this as I only have two areas. From what I have read, it seems that I either have to use an ANCOVA/GLM/mixed-effect model, where I define PA status as both a random effect and a fixed effect. I intended to nest sites within PA, but it seems that as there is only one datapoint per site it will not work as a nested object.

For those familiar with R, here are some codes I have tried:

pcoaPAanovadata1<-read.csv("PCoA\\data\\combined data PCoA axis 1 with distance variables.csv", header=T)

> str(pcoaPAanovadata1)
'data.frame': 80 obs. of 7 variables:
$ PCOA : num -0.2215 -0.3521 -0.0611 0.3434 -0.3624 ...
$ PA.stat: Factor w/ 2 levels "N","P": 1 1 1 1 1 1 1 1 1 1 ...
$ village: num 33.6 33.7 39.9 37.9 34 ...
$ road : num 4.18 3.8 0.89 0.1 3.43 5.49 1.86 5.04 0.79 0.88 ...
$ track : num 8.11 6.48 3.11 2.71 4.49 5.35 1.25 4.03 7.62 6.77 ...
$ site : Factor w/ 80 levels "M1_11","M1_17",..: 1 2 3 4 5 6 7 8 9 10 ...
$ rich : num 3.27 1.79 7.31 0.82 1.79 1.82 2.45 0.82 5.47 2.79 ...

##compare community composition turnover at different PAs

##below specifies a null model where the slope deviates as a result of the random effect
z0 <- lmer(rich ~ 1, random = ~ 1 | pastat/site, data = pcoaPAanovadata1)
summary(z0)
z1 <- lme(rich ~ pastat, random = ~ 1 | pastat/site, data = pcoaPAanovadata1)
summary(z1)
anova(z0,z1)

##impacts of distance variables
zz <- lme(pcoa ~ road, random = ~ 1 | pastat/site, data = pcoaPAanovadata1)
summary(zz)



The errors I get from the lme(linear mixed effect model):
Warning message:
In pt(-abs(tTable[, "t-value"]), tTable[, "DF"]) : NaNs produced

The error I get from the ANOVA:
Warning message:
In anova.lme(z0, z1) :
fitted objects with different fixed effects. REML comparisons are not meaningful.

Firstly, I was hoping to just clarify whether the test I am running is correct. Secondly, it'd be great if someone could tell me what the errors mean. I apologise if my question is poorly phrased, I am relatively new to R and the statistics I am using.

Any help at all would be greatly appreciated.

ted