Tried p + scale_fill_discrete(name = "New Legend Title") but legend title still not changing

I combined two plots of predicted mixed effect model and trying to change the legend title "sR" so that it is friendly to read but I couldn't get it to work when using p + scale_fill_discrete(name = "New Legend Title"). I believe it should be a simple fix but still scratching my head why I can't get it to work even after reading other posts in stackoverflow. Can someone help me please? Thanks.

Below is my R code for you to reproduce the problem:

library(nlme)
library(ggplot2)
library(ggeffects)

#For lowW dataset:

SurfaceCoverage <- c(0.04,0.08,0.1,0.12,0.15,0.2,0.04,0.08,0.1,0.12,0.15,0.2)
TotalSurfaceEnergy <- c(139.31449,105.17776,105.38411,99.27608,92.29064,91.55114,84.44251,78.40453,74.66656,73.33242,72.42429,77.08666)
sample <- c(1,1,1,1,1,1,2,2,2,2,2,2)

lowW <- data.frame(sample,SurfaceCoverage,TotalSurfaceEnergy)

lowW$sample <- sub("^", "Wettable", lowW$sample)
lowW$RelativeHumidity <- "Low relative humidity"; lowW$group <- "Wettable"
lowW$sR <- paste(lowW$sample,lowW$RelativeHumidity)

dflowW <- data.frame(
  "y"=c(lowW$TotalSurfaceEnergy),
  "x"=c(lowW$SurfaceCoverage),
  "b"=c(lowW$sample),
  "sR"=c(lowW$sR)
)

mixed.lme <- lme(y~log(x),random=~1|b,data=dflowW)
pred.mmlowW <- ggpredict(mixed.lme, terms = c("x"))

#For highW dataset: 
  
SurfaceCoverage <- c(0.02,0.04,0.06,0.08,0.1,0.12,0.02,0.04,0.06,0.08,0.1,0.12)
TotalSurfaceEnergy <- c(66.79554,61.46907,57.56855,54.00953,54.28361,55.15855,50.72314,48.55892,47.41811,43.70885,42.13757,40.55924)
sample <- c(1,1,1,1,1,1,2,2,2,2,2,2)

highW <- data.frame(sample,SurfaceCoverage,TotalSurfaceEnergy)

highW$sample <- sub("^", "Wettable", highW$sample)
highW$RelativeHumidity <- "High relative humidity"; highW$group <- "Wettable"
highW$sR <- paste(highW$sample,highW$RelativeHumidity)

dfhighW <- data.frame(
  "y"=c(highW$TotalSurfaceEnergy),
  "x"=c(highW$SurfaceCoverage),
  "b"=c(highW$sample),
  "sR"=c(highW$sR)
)

mixed.lme <- lme(y~log(x),random=~1|b,data=dfhighW)
pred.mmhighW <- ggpredict(mixed.lme, terms = c("x"))

# Combine two predicted mixed effect model into a single ggplot:

p <- ggplot() +
  #lowa plot
  geom_line(data=pred.mmlowW, aes(x = x, y = predicted)) +          # slope
  geom_ribbon(data=pred.mmlowW, aes(x = x, ymin = predicted - std.error, ymax = predicted + std.error), 
              fill = "lightgrey", alpha = 0.5) +  # error band
  geom_point(data = dflowW,                      # adding the raw data (scaled values)
             aes(x = x, y = y, shape = sR)) + 
  
  #higha plot
  geom_line(data=pred.mmhighW, aes(x = x, y = predicted)) +          # slope
  geom_ribbon(data=pred.mmhighW, aes(x = x, ymin = predicted - std.error, ymax = predicted + std.error), 
              fill = "lightgrey", alpha = 0.5) +  # error band
  geom_point(data = dfhighW,                      # adding the raw data (scaled values)
             aes(x = x, y = y, shape = sR)) + 
  xlim(0.01,0.2) + 
  ylim(30,150) +
  labs(title = "") + 
  ylab(bquote('Total Surface Energy ' (mJ/m^2))) +
  xlab(bquote('Surface Coverage ' (n/n[m]) )) +
  theme_minimal()

print(p)

p1 <- p + scale_fill_discrete(name = "New Legend Title")

print(p1)