# How to generate Beta diversity boxplot from phyloseq object?

I am working with profiled metagenomic taxonomy abundance data. I want to generate beta diversity (bray-curtis) boxplot from a phyloseq object where two groups (control and test) will be shown. Something like this figure:

How can I make it? Can anyone share any code or tutorial?

Thanks, dpc

bray_for_plot <- data.frame(phylo_obj@bray)