How is pairwise PERMANOVA/adonis a valid non-parametric approach for pairwise comparisons

Assume that we have taken independent random samples of several individuals from 5 locations that represent 5 populations. The design is fairly unbalanced: the number of individuals sampled from each location is variable (example in code). For each individual, we measured some continuous random variable (e.g., the concentration of some chemical, assuming its value is purely a function of where the respective individual was sampled), and we wish to understand which locations are different using this variable (the concentration of the chemical).
I will simulate example of these data here:

set.seed(123)
data <- data.frame(group = factor(rep(c(paste0("G",1:5)), c(10,24,10,12,9))),
                   val = c(rnorm(10, mean=1.34,sd=0.17), 
                           rnorm(24, mean = 1.14, sd=0.11),
                           rnorm(10, mean=1.19, sd=0.15),
                           rnorm(12, mean=1.06, sd=0.11),
                           rnorm(9, mean=1.09, sd = 0.10)))

Let us assume that group denotes the 5 sampling locations, and val denotes the concentration of the chemical in each sample.

Assume that we have conducted a series of omnibus non-parametric tests (using distance matrices, PERMANOVA) that suggest there are differences between the concentration of the chemical between the groups. We want to conduct pairwise analyses to see which groups are different, so we calculate a dissimilarity matrix (using the vegan package in R for this example) for these data using Euclidean distance:

library(vegan)
dist.mat <- vegdist(data$val, method = "euclidean")

There are several packages out there with functions that conduct pairwise versions of what is generally referred to as PERMANOVA/ permutaional MANOVA/ adonis test (based on the approaches by Marti J. Anderson: Anderson, M.J. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecology, 26: 32-46.)
Some examples can be found in these links:
https://www.rdocumentation.org/packages/RVAideMemoire/versions/0.9-77/topics/pairwise.perm.manova
https://github.com/pmartinezarbizu/pairwiseAdonis/blob/master/pairwiseAdonis/R/pairwise.adonis.R
https://rdrr.io/github/gauravsk/ranacapa/man/pairwise_adonis.html
https://rdrr.io/github/GuillemSalazar/EcolUtils/man/adonis.pair.html

Most of the examples I run across that demonstrate how to use and interpret these approaches in R are applied to count data (mostly species composition in ecological data, which corresponds with Anderson’s intentions). However, there are examples of these approaches being applied in other situations, notably in my interest, for continuous data, such as that presented above. In situations like the one I have proposed, if we calculate a dissimilarity matrix using Euclidean distance, and perform the pairwise version of this procedure (I will use the pairwise.adonis() function by @pmartinezarbizu, 2nd link above):

library(pairwiseAdonis)
#default is 999 permutations
res<-pairwiseAdonis::pairwise.adonis(dis.mat, data[,"group"])
res[,3:5] <- round(res[,3:5],2)  
res

      pairs Df SumsOfSqs F.Model   R2 p.value p.adjusted sig
1  G1 vs G2  1      0.26   19.95 0.38   0.001       0.01   *
2  G1 vs G3  1      0.09    5.06 0.22   0.038       0.38    
3  G1 vs G4  1      0.28   23.98 0.55   0.001       0.01   *
4  G1 vs G5  1      0.34   18.72 0.52   0.001       0.01   *
5  G2 vs G3  1      0.02    1.78 0.05   0.172       1.00    
6  G2 vs G4  1      0.01    0.95 0.03   0.323       1.00    
7  G2 vs G5  1      0.04    2.92 0.09   0.093       0.93    
8  G3 vs G4  1      0.05    3.87 0.16   0.077       0.77    
9  G3 vs G5  1      0.09    4.62 0.21   0.047       0.47    
10 G4 vs G5  1      0.01    0.78 0.04   0.376       1.00    

Am I incorrect in saying that this is equivalent to a simple series of pairwise anovas with p-values calculated according to the observed F statistics probability under the empirical null distribution that was generated through random permutations of group membership (or "location" membership in this case). If so, how can this be a valid non-parametric approach to pairwise comparisons? Let me explain:

Lets use groups 1 and 5 (G1 and G5) as an example:

library(dplyr)
ex <- c("G1","G5")
se <- function(x) sd(x) / sqrt(length(x))
data%>%
  dplyr::filter(., group %in% ex)%>%
  group_by(group)%>%
  summarise_at(., "val", list(mean=mean,med=median,sd=sd,se=se))%>%
  mutate(across(is.numeric, round, 2))

# A tibble: 2 x 5
  group  mean   med    sd    se
  <fct> <dbl> <dbl> <dbl> <dbl>
1 G1     1.35  1.33  0.16  0.05
2 G5     1.05  1.06  0.07  0.02

we know they were heteroscedastic to begin with, but still lets fit a lm and look at residual plots:

ex <- c("G1","G5")
dat2 <- data%>%dplyr::filter(., group %in% ex)
summary(lm(val~group, dat2))
plot(lm(val~group, dat2))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.22775 -0.08020 -0.00070  0.06125  0.27887 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.35269    0.04061  33.306  < 2e-16 ***
groupG5     -0.29792    0.05901  -5.049  9.9e-05 ***

Now lets do regular ANOVA:

ex <- c("G1","G5")
dat2 <- data%>%dplyr::filter(., group %in% ex)
summary(aov(val ~ group, data = dat2))

            Df Sum Sq Mean Sq F value  Pr(>F)    
group        1 0.4204  0.4204   25.49 9.9e-05 ***
Residuals   17 0.2804  0.0165    

Now lets do the pairwise adonis procedure for just these two groups:

set.seed(123)
dis.mat <- vegdist(dat2$val, method="euclidean")
res<-pairwiseAdonis::pairwise.adonis(dis.mat, dat2[,"group"])
res[,3:5] <- round(res[,3:5],2)  
res

     pairs Df SumsOfSqs F.Model  R2 p.value p.adjusted sig
1 G1 vs G5  1      0.42   25.49 0.6   0.001      0.001  **

We get the same observed model, as we should, and a p-value that is based off of permutations. A little look under the hood shows me that all this function really does is perform the adonis function on each pair (or the single pair in this case). To demonstrate, we get the same answer by doing this:

set.seed(123)
dis.mat <- vegdist(dat2$val, method="euclidean")
res2<- adonis(dis.mat ~ dat2$group, method = "euclidean")
res2$aov.tab

Terms added sequentially (first to last)

           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
dat2$group  1   0.42041 0.42041  25.488 0.59989  0.001 ***
Residuals  17   0.28041 0.01649         0.40011           
Total      18   0.70082                 1.00000     

So all we have really done here (still using the G1vsG5 example) is

1. calculate an F test for the raw comparison (of the dissimilarity values, which are mathematically equivalent to differences in the raw data since we used Euclidean distance),
2. shuffled the raw data and calculated a new F statistic
3. repeated step 2 999 times to create an empirical F distribution (generating the null model),
4. and finally calculated the probability of the observed F value occurring under the null model:

perm<-permustats(res2)
densityplot(perm)

So if all we are really doing is comparing F values that are calculated from the normal anova/linear model, is this really a valid "non-parametric" approach to make pairwise comparisons?