I have a data set where community composition data was collected from 2 sites in a paired design: under the canopy and nearby canopy. We used mvabund package in R to check whether the communities differed. They do differ. The main problem is to find a way to plot the data and show that communities are different, accounting for the paired design. We plotted the data through nMDS, but the communities did not seem to be different at all. I'm trying to use boral package but I haven't found in any example how to control for the pairing (blocks).
The updated version of BORAL allows for inclusion of site effects, which you can group as needed (replicates, nested, blocked, etc). See examples from https://www.rdocumentation.org/packages/boral/versions/1.7/topics/boral.
You can use results from BORAL in ggplot to create a figure similar to NMDS (using latent variables as axes), to show community differences.