Look at the chunk of ggplot code underneath the first ordination. I’m stuck at the fit env variables as vectors part. I don’t know what to do to get veg.sites and veg.env to be the same length. In your example, you’re using dune, dune.env, and dune.sites, which are all the same length.

## Warning: package 'bindrcpp' was built under R version 3.4.3
## Call: capscale(formula = veg.d[3:14] ~ 1, distance = "bray")
## 
##               Inertia Rank
## Total          17.950     
## Unconstrained  25.491   71
## Imaginary      -7.541  166
## Inertia is squared Bray distance 
## Species scores projected from '[' 'veg.d' '3:14' 
## 
## Eigenvalues for unconstrained axes:
##  MDS1  MDS2  MDS3  MDS4  MDS5  MDS6  MDS7  MDS8 
## 6.511 4.094 2.539 1.488 1.127 0.882 0.827 0.819 
## (Showed only 8 of all 71 unconstrained eigenvalues)
## Warning in match.fun(FUN)(...): "border" is not a graphical parameter
## Warning in match.fun(FUN)(...): "border" is not a graphical parameter

## Warning in match.fun(FUN)(...): "border" is not a graphical parameter

## Warning in match.fun(FUN)(...): "border" is not a graphical parameter
## Warning in rep(fill, length = nrow(x)): 'x' is NULL so the result will be
## NULL

pacman::p_load(tidyverse, plyr, dplyr, vegan)


veg.cap <- capscale(veg.d[3:14]~1,  data = veg.d)

# Extract scores from ord results
veg.spp <- scores(veg.cap, scaling = 1)$species %>%
  as.data.frame() %>%
  rownames_to_column("species")

veg.sites <- scores(veg.cap, scaling = 2)$sites %>%
  as.data.frame() 

veg.env <-
  veg.all %>%
  
  mutate(management=case_when(
    Location=="REF"~"continuous",
    Location=="WAG"~"twice-over",
    Location %in% c("BOB","BAR") &
      year==2018 & Year==2018 & Season=="Spring"~"burned",
    year==2017 & Year==2017 & Season=="Spring"~"burned",
    TRUE~"unburned")) %>%
  
  select(year, management, BareGrd:Litter_4)



# fit env variables: vectors
veg.vec <- scores(envfit(cbind(veg.sites$MDS1, 
                                veg.sites$MDS2),  
                          veg.env,  
                          choices=c(1:2)), 
                   "vectors")  %>%
  round(., 3)  %>%
  as.data.frame() 

nrow(veg.sites) #238
nrow(veg.env) #4436
nrow(veg.all) #4436
nrow(veg.d) #238