df <- mtcars

fit1 <- gls(mpg ~ cyl + disp, data = df)
fit2 <- gls(mpg ~ cyl + disp + hp, data = df)

fit1ml <- update(fit1, . ~ ., method = "ML")
fit2ml <- update(fit2, . ~ ., method = "ML")
anova(fit1ml, fit2ml)

> anova(fit1ml, fit2ml)

     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fit1ml     1  4 167.1456 173.0086 -79.57282
fit2ml     2  5 168.0184 175.3471 -79.00921 1 vs 2 1.127214  0.2884