While I feel pretty confident that you can use the predictions from nlme::gls() as the predictions of the response and not the linearly transformed predictions that was done for the GLS algorithm, I think finding that info will take some sleuthing.
The Pinheiro and Bates book, Mixed Effects Models in S and S-plus, is likely a good place to start looking fore more info. This is the book that package authors wrote about nlme, often considered the "go-to" for understanding the package.
Another idea is to look through the gls() function code. While reading code can be difficult (for me, at least), if you can wade through it you may well be able to see if results (betas, fitted values) are on the original or linearly transformed scale. See the function code by running, e.g., gls() in your R Console. Things look fairly complicated to me so this might not go to far.
My final idea is to create data where you know what the predictions should be on the original scale and then see if gls() recreates them. This could even be something simpler like fitting a model with OLS via lm() and then gls() with heteroskedastic errors argument and see if the predictions change.