I am trying to make comparisons between treatments of my heritability estimates, values ​​that I calculated with the following function:

data<- as.data.frame(treatment=c("control","ivermectin","pla2guicide","parasitism"),temperature=c("29","33"), pronotum=c(4.12,4.21,4.33,4.26,4.34,4.21,4.52,4.35))
data

her_pro <- function(data, indices) {
mod <- lmer(pronotum ~ 1 + (1|terrarium), data = data[indices, ])
df <- as.data.frame(VarCorr(mod))
family<-df[1,4]
residual<-df[2, 4]
heritability <- family / (family + residual)
return(heritability)
}

#I used that function to obtain a heritability index, from a mixed model.
#Then I made a Bootstrap as follows:

BC29_p <- boot(data=data, statistic=her_pro, R=100)

Now, I want to know if that value of heritability is statistically different between treatments.

What test can I do to compare heritability values?
I’ve been trying with the finalfit library and the boot_predict function, but I’m not sure if I’m doing it right?