Why does the order of factors (order of samples in dataset) affect PERMANOVA results with adonis2?

Question

I am trying to run PERMANOVA tests on multiomics datasets collected from coral samples using adonis2 (v 2.6-5). I found when the order of samples in my dataset changes, the P-values change, sometimes markedly. For instance, I have 4 time points (0, 1, 3, 5), and if the order of samples is 0, 1, 3, 5 the PERMANOVA results are different than if the samples are ordered TP 1, 3, 5, 0.

I have made sure the order of the samples in the sample information with the factors (master dataframe) matches the sample order in the data frame, both below for reference. The data frame has the samples as rows and transcript names as columns.

When the samples are in this order the PERMANOVA results are as follows.

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  strata 
Permutation: free
Number of permutations: 999

adonis2(formula = data.dist ~ TP * Treatment, data = fac.data, permutations = 999, method = "bray", strata = fac.data$Tank)
             Df SumOfSqs      R2      F Pr(>F)  
TP            3 0.040207 0.28861 2.8343  0.057 .
Treatment     1 0.020644 0.14819 4.3658  0.101  
TP:Treatment  2 0.012260 0.08800 1.2963  0.245  
Residual     14 0.066200 0.47520                
Total        20 0.139311 1.00000

The PERMANOVA results for when the samples are reordered (in both data frames) to TP 1, 3, 5, 0 (but nothing else changes) are

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  strata 
Permutation: free
Number of permutations: 999

adonis2(formula = data.dist ~ TP * Treatment, data = master, permutations = 999, method = "bray", strata = fac.data$Tank)
             Df SumOfSqs      R2      F Pr(>F)   
TP            3 0.040207 0.28861 2.8343  0.002 **
Treatment     1 0.020644 0.14819 4.3658  0.010 **
TP:Treatment  2 0.012260 0.08800 1.2963  0.231   
Residual     14 0.066200 0.47520                 
Total        20 0.139311 1.00000

I will attach my relevant code below, but in summary, for TP 1, 3, and 5 samples were collected under high temperature and ambient temperatures (control) and grown in tanks. TP 0 are field samples collected from the same coral colonies to get a baseline of any tank effects. I am analyzing proteomic, metabolomic, transcriptomic, and 16S amplicon sequence data and treating each the same after their respective normalization. I am additionally taking the square root of each value and then creating a dissimilarity matrix.

This is only the code pertaining to the PERMANOVA test, which is part of a larger analysis, which is why there are so many packages, but I'll show each just in case it matters. The R version is 4.2.0.

# Load libraries
library(magrittr)
library(lubridate)
library(tidyverse)
library(seacarb) 
library(matrixStats)
library(vegan)
library(RVAideMemoire)
library(lme4)
library(ape)
library(emmeans)
library(gridExtra)
library(multcompView)
library(plotrix)
library(reshape2)
library(ggpubr)
library(sva, warn.conflicts = FALSE)
library(mixOmics)
library(ggforce)
library(dplyr)
library(plotly)
library(gapminder)
library(htmlwidgets)

# Set seed
set.seed(54321)

### Format data for downstream analyses ###

# Filter for data with total counts >= ____
filtered.data <- data.df %>% 
            replace(is.na(.), 0) %>%
            dplyr::mutate(sum = rowSums(across(where(is.numeric))))
filtered.data <- filtered.data[!(filtered.data$sum <= 100),]
filtered.data$sum <- NULL

### Format samplesInfo for downstream analyses ###
master <- master %>%
            mutate_all(as.character)
master$TP <- as.factor(master$TP)
master$Treatment <- as.factor(master$Treatment)
master$Tank <- as.factor(master$Tank)

str(master)


### PERMANOVA ###

# Transpose df
t.filtered.data <- t(filtered.data)
t.filtered.data <- as.data.frame(t.filtered.data)

# Complete dataframe for PERMANOVA
fac.data <- merge(master, t.filtered.data, by = 0)
rownames(fac.data) <- fac.data[[1]]
fac.data[[1]] <- NULL

# Use square root or proportions to minimize influence of most abundant groups
sum(is.na(t.filtered.data))
data.mat <- sqrt(data.mat)
sum(is.na(data.mat))
which(is.na(data.mat), arr.ind = TRUE) #where NA's are located, if present

# Create a dissimilarity matrix (vegan)
data.dist <- vegdist(data.mat, method = 'bray')

# Run perMANOVA (vegan)
mod.data <- adonis2(data.dist ~  Treatment * TP, data = fac.data, permutations = 999, strata = fac.data$Tank, method = 'bray')
mod.data

If anyone could please help me understand what is happening and how to correct for this so I get accurate PERMANOVA results, I would appreciate it greatly! Also, if I haven't included any pertinent information, please let me know.

Answer 1

Please read the documentation: it is intended to answer potential questions. In this case, you should go to argument by. Quoting from help(adonis2):

by: ‘by = "terms"’ will assess significance for each term (sequentially from first to last), setting ‘by = "margin"’ will assess the marginal effects of the terms (each marginal term analysed in a model with all other variables), and ‘by = NULL’ will assess the overall significance of all terms together. The arguments is passed on to ‘anova.cca’.

The default is to analyse terms sequentially which is the same as in the standard anova of base R.