Problems using glmulti package in RStudio to predict best model

Question

So I am trying to use the "glmulti" package find the best combination of variables (or best model) for my response variable Bio_class. Bio_class is a categorical variable that contains 5 outcomes (4 species and 1 absence option). I have 56 numerical variables (temperature, salinity, chlorophyll, etc.) that I am trying to test in the model. I have followed other posts online for how to use the package but I am getting errors. I am fairly new to R so any recommendations on whether to include random effects or not would be appreciated. Please let me know if you have any questions. I know this is ALOT of data to put into the package but I hope it works. Ive also with less variables and excluded substrate and I get a different error.

Many thanks

library(glmulti)
library(metafor)

data.na<- na.omit(alldata) #remove any rows with NAs in them
anyNA(data.na) #FALSE = no NAs in dataset

rma.glmulti <- function(formula, data, ...)
rma(formula, vi, data=data, method="ML", ...) #defining the function 'rma.glmulti'

res <- glmulti(
    Bio_class~Temp_winter+Temp_spring+Temp_summer+Temp_autumn+Temp_yr_avg
    +Temp_yr_max+Temp_yr_min+Sal_winter+Sal_spring+Sal_summer+Sal_autumn
    +Sal_yr_avg+Sal_yr_max+Sal_yr_min+Chl_winter+Chl_spring+Chl_summer+Chl_autumn
    +Chl_yr_avg+Chl_yr_max+Chl_yr_min+Nit_winter+Nit_spring+Nit_summer+Nit_autumn
    +Nit_yr_avg+Nit_yr_max+Nit_yr_min+ZSD_winter+ZSD_spring+ZSD_summer+ZSD_autumn
    +ZSD_yr_avg+ZSD_yr_max+ZSD_yr_min+Phos_winter+Phos_spring+Phos_summer          
    +Phos_autumn+Phos_yr_avg+Phos_yr_max+Phos_yr_min+SWH_winter+SWH_spring
    +SWH_summer+SWH_autumn+SWH_yr_avg+SWH_yr_max+SWH_yr_min+Depth
    +Dist_to_major_cities+Dist_to_complete_cities+Dist_to_port
    +Dist_to_major_river+Dist_to_complete_river+Dist_to_coast, 
    data=data.na,
    level=1,
    fitfunction=rma.glmulti, 
    crit="aicc", 
    confsetsize=7.2057594e+16
)

print(res)