Thanks to @dipetkov's hints in the comments, I was able to correct my mistakes, and reproduce the calculation as follows:

# Lower band calculation
xh <- lrmodel$model$wt[1]               # 2.62
n <- length(lrmodel$residuals)          # 32
tval <- qt(p = 0.01/2, df = n - 2)      # -2.749996
xbar <- mean(lrmodel$model$wt)          # 3.21725
mse <- sum(lrmodel$residuals^2)/(n - 2) # 9.277398
yhat <- preds[1, "fit"]                 # 23.28261

# upper band (31.83807)
yhat - tval * sqrt(mse * (1 + 1/n + (xh - xbar)^2/sum((lrmodel$model$wt - xbar)^2)))

# lower band (14.72715)
yhat + tval * sqrt(mse * (1 + 1/n + (xh - xbar)^2/sum((lrmodel$model$wt - xbar)^2)))
```

Thanks to @dipetkov's hints in the comments, I was able to correct my mistakes, and reproduce the calculation as follows:

# Lower band calculation
xh <- lrmodel$model$wt[1]               # 2.62
n <- length(lrmodel$residuals)          # 32
tval <- qt(p = 0.01/2, df = n - 2)      # -2.749996
xbar <- mean(lrmodel$model$wt)          # 3.21725
mse <- sum(lrmodel$residuals^2)/(n - 2) # 9.277398
yhat <- preds[1, "fit"]                 # 23.28261

# upper band (31.83807)
yhat - tval * sqrt(mse * (1 + 1/n + (xh - xbar)^2/sum((lrmodel$model$wt - xbar)^2)))

# lower band (14.72715)
yhat + tval * sqrt(mse * (1 + 1/n + (xh - xbar)^2/sum((lrmodel$model$wt - xbar)^2)))