How can i get predictions with CI from lmerTest models?

We are currently working with plant phenology.

We built a linear mixed model for each species present in the study area.

We set Days From Snowmelt (The sum of days from snowmelt to the visit day along the summer) as the response variable while Mean phenology (mean phenology state for each plot ( there are 3 on each locality) is calculated by the mean phenological state from the 12 subplots into each plot is divided. from 1-6, the higher the number the more advanced the cycle). year and plot nested within the locality are set as random factors.

Once the model is built and revised, we want to predict the days from snowmelt for each species to achieve the phenological phases of interest, which happen to have a mean of 2, 3, 4, and 5. (corresponding to vegetative, flowering, fruit development and dispersion, respectively) I have tried the function predict() but I get no heterogeneity between phases for each species, the progression seems to be linear (as shown in the image file).

Could this be just because is a linear model so will it only give linear responses? Are there any other ways to get predictions from these kinds of models and show their CI?

How can i get predictions with CI from lmerTest models?

I think you probably mean pediction intervals. You can use the predictInterval function in the merTools package. For example:

library(lmerTest); library(merTools)

fm1 <- lmer(Reaction ~ Days + (Days|Subject), data = sleepstudy)
head(predictInterval(fm1, level = 0.95, seed = 123, n.sims = 100))

Could this be just because is a linear model so will it only give linear responses?

Yes ! If you fit a linear model, then the predictions will be linear. Of course, you can model nonlinearity with a linear model in several ways including transformation(s), nonlinear terms (the model is still linear in the parameters) and splines.