How to [properly] correlate environmental variables to NMDS ordination axes?

Question

Data + setup:

• I've constructed an NMDS ordination from a Bray-Curtis distance matrix calculated from relativized abundances (basal areas) of trees. Samples include ~40 forested plots that have been resampled 13-15 times each.

• I also have ~60 environmental variables (soil, topography, etc.) that were measured from each plot in the last sampling year. Many variables relate and/or are correlated to one another.

Goal:

• I want to determine if any of my environmental variables explain the NMDS structure and if so to what degree.

Attempt:

• I opted to use vf() from the ecodist package in R (though I know a lot of people use envfit() or ordisurf functions from the vegan package). - From the help page, vf provides the following:

  matrix with the first 2 columns containing the scores for every variable in each of the 2 dimensions of the ordination space. r is the maximum correlation of the variable with the ordination space, and pval is the result of the permutation test.

• Actually, all of these functions generate either a Pearson's r or an R².

Concern

• Many of my environmental variables are non-normally distributed (and are not made normal via log transformation).

  • As tested with QQ plots, histograms, and Shapiro-Wilk test

• My NMDS ordination scores (i.e., the sample point coordinates in my NMDS plot) are themselves not normally distributed (but of course, why would they be?!)

QUESTIONS:

1. Is a Pearson correlation appropriate for assessing the correlation between any [non-normal] environmental variable and the NMDS axes?

• If so, do I have to transform any data to make it work??

  • If yes, I don't think it makes sense to transform the ordination axes, so what do I do?

2. Given #1, What is the appropriate way to determine these correlations?