Why is constrained RDA symmetrical? I am having a very weird behaviour when performing constrained RDA. The hypothesis I am testing is: is there a difference between treated and control samples? The experimental design comprise 10 sites, where a control and a treatment were sampled.
As far as I understand, this is the classical example of random effects caused by the sampling site so my model would be:
data_agglomerated_at_taxa_level ~ Treatment + Condition(Site)
However, something very interesting happens: as you can see from the picture:

what I am observing are samples that are symmetrical (I and III quadrants, II and IV quadrants). Each couple of symmetrical data point is the couple $(control_{s}, treatment_{s})$, $s = site$.
If I remove the constraint on the Site, then the data points are not symmetrical anymore so, as far as I can tell, there is something happening with the Site. So said, I can't really understand what that might be. Any ideas?
PS: i am performing the analysis with Phyloseq ordinate() function and the issue seems to be related to the type of data in the sample_data(phylo_object).
UPDATE:
i am adding a reproducible example, and another issue that i found:
library("vegan")

# create fake otu
otu_fake <- matrix(0, nrow=10, ncol=100)

# zero inflated random fake otu data
for (i in 1:100) {
    zero_infl <- ifelse(rbinom(n = 10, size = 1, prob = 0.90)==1, 0, 1)
    ones <- which(zero_infl==1)
    otu_fake[ones, i] <- sample.int(10, length(ones))
}
colnames(otu_fake) <- paste0("otu_", 1:100)
rownames(otu_fake) <- paste0("sample_", 1:10)

# create fake meta
fake_meta <- data.frame(sample=paste0("sample_", 1:10), treatment=c(rep("CTR", 5), rep("TRT", 5)), site=c(1:5, 1:5))

# perform RDA
rda_veg <- rda(formula = otu_fake ~ treatment+Condition(site), data = fake_meta)

# plot it
dev.new()
plot(rda_veg, main="Site as integer")

# cast site to factors
fake_meta$site <- as.factor(fake_meta$site)

# perform RDA
rda_veg <- rda(formula = otu_fake ~ treatment+Condition(site), data = fake_meta)

# plot it
dev.new()
plot(rda_veg, main="Site as factors")

This is an example of the experiment i am working on. The funny thing is that the symmetry appears only after casting the site variable to factor (or character).
So, updated question is: is this due to the design/mathematical structure of the data, is there something else, or am I missing something?
UPDATE on the update
I provide here more code where the symmetry issue can be tested, while the factor/integer problem remains.
To see the symmetry, samples needs to equal to 10. For higher values, the symmetry is lost.
library("vegan")

# choose a number of taxa and samples
taxas <- 100
samples <- 30

# with samples = 10 the setting is:
# 1 site comprise 1 control and 1 treated sample

# with samples = 20, 30, 40... the setting becomes:
# 1 site comprise 2 controls and 2 treated samples

# create fake otu
otu_fake <- matrix(0, nrow=samples, ncol=taxas)

# zero inflated random fake otu data
for (i in 1:taxas) {
    zero_infl <- ifelse(rbinom(n = samples, size = 1, prob = 0.90)==1, 0, 1)
    ones <- which(zero_infl==1)
    otu_fake[ones, i] <- sample.int(samples, length(ones))
}
colnames(otu_fake) <- paste0("otu_", 1:taxas)
rownames(otu_fake) <- paste0("sample_", 1:samples)

# create fake meta
fake_meta <- data.frame(sample=paste0("sample_", 1:samples), treatment=c(rep("CTR", samples/2), rep("TRT", samples/2)), site=1:5)

# perform RDA
rda_veg <- rda(formula = otu_fake ~ treatment+Condition(site), data = fake_meta)

# plot it
dev.new()
plot(rda_veg, main=paste("Site as integers with ", samples, " samples"))

# cast site to factors
fake_meta$site <- as.factor(fake_meta$site)

# perform RDA
rda_veg <- rda(formula = otu_fake ~ treatment+Condition(site), data = fake_meta)

# plot it
dev.new()
plot(rda_veg, main=paste("Site as factors with ", samples, " samples"))


I don't know much about RDA and CCA, but this might help you into the right direction.
The function seems to apply some transformation before performing PCA. This transformation seems to center the data in a way such that the average PC score per site is zero.
### rda uses Ybar, some centered and re-scaled values per column
### the object contains 'rda_veg[6]$pCCA$Fit' this is a vector
### that resembles the center of the pca scores per site
### when you subtract it then you get the symmetry

### perform pca with this
otu_transform = rda_veg[3]$Ybar-rda_veg[6]$pCCA$Fit


pc = prcomp(otu_transform, center = FALSE)$rotation

### plot pc1 and pc2
### still not exactly the same but similar
sites1 = otu_transform[,] %*% pc[,1]
sites2 = otu_transform[,] %*% pc[,2] 
plot(sites1,sites2, col = fake_meta$site)


The upper graph pc1 vs pc2 is made with your code but adding a choices parameter plot(rda_veg, choices = c(2,3), main="Site as factors") 
