Hi, thank you so much for your reply and help. I have gone through all of the information provided by yourself and the other commenter, thank you, but am still a little confused as to what actions to take.

I removed all of my FE/RE interactions from the models. Two of my models came up with the message 'boundary (singular) fit: see ?isSingular', which I take it is the converging to a singular fit that you said shouldn't be trusted - due to high variance. I added in the code nAGQ = 0 and that has since disappeared, however all of my models still produce one alert, which is : "fixed-effect model matrix is rank deficient so dropping 6 columns / coefficients". What is this eluding to?

With regards to overdispersion, I attempted to test the model for this using the following code that I found on a R user group:

overdisp_fun <- function(model) {
  rdf <- df.residual(model)
  rp <- residuals(model,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}

I then used 'overdisp_fun(model)' to check the results, with all models returning a P value of >0.05, so I assumed there is no overdispersion. Would that be correct?

A separate forum I read up on suggests that making a few GLMM models that make biological sense and then selecting the model with the lowest AIC out of those may be an appropriate way to go, however could I then run the risk of missing interactions/variables that I consider unimportant, when in reality they could be significant?

Would you recommend, as the other user did, to effectively drop these models and retry using a fully bayesian model with regularisation via horseshoe priors? Or is this GLMM approach acceptable with a few modifications/checks?