corr.axes output

Hi,
I have been running corr.axes today using a a metadata input file with chemistry data, and after plotting my data in R, I am not sure if my output is correct, or that I am missing something.
The command line was the following:

orr.axes(axes=final.an.thetayc.0.03.lt.ave.nmds.axes, metadata=OFS_metadata.txt, method=spearman, numaxes=2, label=0.03)

My chemistry data looks like this (I reduced it to make my point):

group TC TOC TN
PM10A04 3.66 2.53 0.209
PM10B04 3.78 2.56 0.308
PM10C04 3.76 2.54 0.327
PM11A04 3.65 2.46 0.331
PM11B04 2.98 2.03 0.185
PM11C04 3.8 2.43 0.122
PM12A04 3.53 2.36 0.152
PM12B04 3.54 2.38 0.196
PM12C04 3.45 2.4 0.170
RDA04 3.45 2.4 0.139
RDB04 3.73 2.28 0.117
RDC04 3.4 2.35 0.111
REA04 3.3 2.31 0.177
REB04 3.59 2.53 0.234
REC04 3.77 2.54 0.120
PM10A40 1.29 0.84 0.654
PM10B40 1.24 0.82 0.776
PM10C40 1.48 1.03 0.659
PM11A40 1.23 0.81 0.713
PM11B40 1.07 0.55 0.482
PM11C40 1.07 0.58 0.518
PM12A40 1.36 0.83 0.501
PM12B40 1.23 0.84 0.564
PM12C40 1.24 0.74 0.556
RDA40 1.63 1.19 0.503
RDB40 1.45 1.07 0.591
RDC40 1.41 1.09 0.602
REA40 1.51 1.1 0.557
REB40 1.55 1.16 0.466
REC40 1.58 1.15 0.623

So just eyeballing the data I see that the TC/TOC and TN data have a different direction with the top half of the data set having higher values for TC/TOC and lower for TN and vice versa.
After running the corr.axes command I get the following output in the file: OFS_metadata.spearman.corr.axes

Feature axis1 p-value axis2 p-value length
TOC -0.754507 0.000001 0.702203 0.000015 1.030714
TC -0.727576 0.000005 0.709771 0.000011 1.016436
TN 0.486096 0.00646 -0.714794 0.000009 0.864419

Plotting this over the NMDS plot using the script found here: http://mothur.ltcmp.net/t/corr-axes-length-meaning-and-r-script-for-biplot/598/3
I get arrows that seem to point in a different direction than I expected. TN is correlated in my plot with the surface samples, which should not be, since the higher TN concentrations are in the deeper samples. So is the output of corr.axes correct, or am I missing something essential here?

Here is the data for my nmds plot, which I edited to sort it:

sample axis1 axis2 depth location
PM10A04 0.294035 -0.331402 4 pockmark
PM10B04 0.279339 -0.349651 4 pockmark
PM10C04 0.223797 -0.322045 4 pockmark
PM11A04 0.208262 -0.363415 4 pockmark
PM11B04 0.32964 -0.352292 4 pockmark
PM11C04 0.252192 -0.388492 4 pockmark
PM12A04 0.158514 -0.356438 4 pockmark
PM12B04 0.051985 -0.417124 4 pockmark
PM12C04 0.149395 -0.287356 4 pockmark
RDA04 -0.005351 -0.333037 4 normal
RDB04 -0.048915 -0.314942 4 normal
RDC04 -0.001089 -0.366386 4 normal
REA04 0.085492 -0.313089 4 normal
REB04 0.060524 -0.376752 4 normal
REC04 0.106204 -0.357482 4 normal
PM10A40 -0.347978 0.334938 40 pockmark
PM10B40 -0.130239 0.119286 40 pockmark
PM10C40 -0.240211 0.335685 40 pockmark
PM11A40 -0.408563 0.288608 40 pockmark
PM11B40 -0.446157 -0.058342 40 pockmark
PM11C40 -0.232386 0.451806 40 pockmark
PM12A40 -0.105009 0.351927 40 pockmark
PM12B40 -0.521901 0.269354 40 pockmark
PM12C40 -0.446411 0.410278 40 pockmark
RDA40 0.16082 0.383536 40 normal
RDB40 0.051326 0.437175 40 normal
RDC40 0.090821 0.43928 40 normal
REA40 -0.006591 0.49465 40 normal
REB40 0.319814 0.491846 40 normal
REC40 0.118639 0.479876 40 normal

Hey coming back to my previous posting.

I calculated in excel the Spearman rank correlation for the ndms axis of my dataset together with the TN factor, using the formula found here:

I rank the data from smallest to the largest value.

If I then calculate the correlation coefficients for TN I get the following.

axis 1 axis 2
TN -0.45895 0.71657

These values are the opposite of what I found yesterday coming out of corr.axes. The values are not identical, but since they are in the same range as the values in corr.axes, I think I can trust my calculations.
When I invert the ranking, with the largest having the smallest rank, I get correlation coefficients with opposite directions and similar as corr.axes, which seems not correct to me.

Is this the right conclusion I draw here?

What version of mothur are you using?

When I run the samples you provided above in 1.29.1, I am getting:

Feature axis1 p-value axis2 p-value length
TC 0.712887 0.000010 -0.700423 0.000016 0.999400
TOC 0.738927 0.000003 -0.694414 0.000021 1.014014
TN -0.458954 0.010739 0.716574 0.000008 0.850951