From reviewer:
A much simpler and standard way to do this is to just put the technical covariates directly into a sequential (Type I) PERMANOVA model of the original, unprojected Wasserstein distance matrix:
adonis2(wdist ~ log(ncod) + log(treelen) + family + exposed * enveloped, data=mdat, by='terms')
This way, the model attributes variance to log(ncod) and log(treelen) first, controlling for them before testing the biological factors, all while using the full distance matrix.
When I run this sequential model on the 205 valid proteins, log(ncod) explains 58.3% of the variance, log(treelen) explains 2.8%, and family explains 12.9%. The biological factors explain a very small slice: exposure explains 0.36% (marginally significant, P = 0.082) and the interaction term explains 0.56% (statistically significant, P = 0.025).
From reviewer: