Use vegdist function before metaMDS? I am trying to do calculate beta-diversity of bacterial community data. I have been following numerous tutorials regarding this subject and I am getting mixed advice. I am trying to calculate beta-diversity and illustrate using NMDS plot.
Some people would say that you have to calculate beta-diversity first using the vegdist function from the vegan package like this:
bray_dist <- vegdist(OTU, method = "bray")
beta_nmds_bray = metaMDS(bray_dist, distance="bray", k=2)

However, there are people saying that vegdist is included in the metaMDS function and it can be calculated like this:
beta_nmds_bray = metaMDS(OTU, distance="bray", k=2)

Naturally I am getting different results and am confused about which route is correct.
I also read that this approach is also representitive of beta-diversity, but not sure. I think it is part of PERMADISP?
bray_dist <- vegdist(OTU, method = "bray")
disp.age_bray <- betadisper(bray_dist, meta_data$stream)
anova(disp.age_bray)
permutest(disp.age_bray, pairwise=TRUE, permutations=1000)
plot(disp.age_bray)

Any advice would be appreciated!
 A: As NMDS is prone to converging to locally optimal solutions, not globally optimal ones, the metaMDS() helpers fit the NMDS from multiple random starts and then compare the fits of subsequent runs. If you get two runs that converge to essentially the same solution (assessed via procrustes rotation) then the metaMDS wrapper will declare convergence and terminate the random starts. If you don't it keeps running until you hit the pre-defined number of fits (set via argument trymax).
The key word above is random; unless you fix the random seed before repeated calls to metaMDS() you'll get different random starts and potentially convergence to different results. That should worry you unless the differences are slight as it implies different local optima that are relatively broad so they keep being found.
This isn't the only thing that metaMDS() does however and it is these extra things that the wrapper does when provided with community data that it doesn't do when provided with a dissimilarity matrix that is causing additional differences.
library('vegan')
data(varespec)

dij <- vegdist(varespec, method = 'bray')

set.seed(1)
ord1 <- metaMDS(varespec, distance = 'bray', k = 2)
set.seed(1)
ord2 <- metaMDS(dij, k = 2)
set.seed(1)
ord3 <- metaMDS(varespec, distance = 'bray', k = 2, autotransform = FALSE, 
                halfchange = FALSE)

If you run that you'll see that the scores for ord1 and ord2 are different
> all.equal(scores(ord1, 'sites'), scores(ord2, 'sites'))
[1] "Mean relative difference: 0.5419395"

even though we fixed the seed. This is the result of the extra help things that metaMDS does. If we turn this off (autotransform = FALSE) and halfchange off and compare ord2 and ord3 then we get the same site scores
> all.equal(scores(ord2, 'sites'), scores(ord3, 'sites'))
[1] TRUE

Although it is worth noting that the objects ord2 and ord3 are themselves not the same — there are other things that we can do if given the community data that we can't do (or don't do) if only given the dissimilarities, but these things don't affect the NMDS solution itself.
With something as complex as metaMDS it really is helpful to read and understand the documentation; not everyone that uses this stuff and writes about it online is aware of the details, which do matter. I'm one of the vegan developers and I had to look some of this stuff up as it's not something I use often in my day job.
Finally, note that computing the Bray-Curtis distances between your points is not the same as computing the beta diversity of your samples. For that to be true you'd have to be equating local and regional diversity as each sample being a region itself.
A: It would really help if you posted some example data.
It shouldn't make any difference whether the input for metaMDS is a dissimilarity matrix you've pre-calculated (bray_dist), or your raw OTU table (OTU). Are you sure you're really getting different results?  When you call beta_nmds_bray do you not see the same 'stress' score?
You will see a slight difference if you try plotting the outputs (beta_nmds_bray) - this is just because the species scores are missing if you use a dissimilarity matrix as your input to metaMDS. If you plot just the 'site' points, the plots should like very similar - identical except for perhaps being on a different scale.
You can calculate a range of beta diversity measures using the betadivers() function, and the betadispers() function can be used to analyse beta diversities with respect to classes or factors.
I recommend reading in full the vignettes for the vegan package:
https://cran.r-project.org/web/packages/vegan/vignettes/diversity-vegan.pdf
https://cran.r-project.org/web/packages/vegan/vignettes/intro-vegan.pdf
