m1 <- lmer(Z_TP~Z_Length + (1|Lake) + (1|Species), data = data, REML=TRUE)
m2 <- lmer(Z_TP~Z_Length + (1+Z_Length|Lake) + (1+Z_Length|Species), data = data, REML=TRUE)
m3 <- lmer(Z_TP~Z_Length + (1|Species), data = data, REML=TRUE)
m4 <- lmer(Z_TP~Z_Length + (1|Lake), data = data, REML=TRUE)
m5 <- lmer(Z_TP~Z_Length + (1+Z_Length|Species), data = data, REML=TRUE)
m6 <- lmer(Z_TP~Z_Length + (1+Z_Length|Lake), data = data, REML=TRUE)
m7 <- lmer(Z_TP~Z_Length + (1+Z_Length|Lake) + (1|Species), data = data, REML=TRUE)
m8 <- lmer(Z_TP~Z_Length + (1|Lake) + (1+Z_Length|Species), data = data, REML=TRUE)

#Bonus model!
M0<-lm(Z_TP~Z_Length,data=data)
# It is always useful to build the simple linear model without any varying intercept and slope factors to see the difference in AICc values. Although "lm" does not use the same estimation method