What is the difference between generalized linear models and generalized least squares Generalized linear models (GLMs) allow the response variable to have arbitrary distributions using a link function
Generalized least squares account for errors that are heteroscedastic and correlated. 
But is generalized least squares a special case of GLMs? or are they two different methods altogether?
 A: No, these are two different things. GLMs are models whose most distinctive characteristic is that it is not the mean of the response but a function of the mean that is made linearly dependent of the predictors. 
GLS is a method of estimation which accounts for structure in the error term. An ordinary linear model could be estimated by GLS if you think that errors are not independent and homoscedastic.
A: As said, those are generally two different things. 
GLM framework = link function + other distribution than normal
GLS framework = generalizes the iid normal in LM to a multivariate normal, which allows specifying correlations between the residuals + change of dispersion (in R, this is easiest done with nlme, which provides several corClasses to specify correlation structures, and the varFun function to specify a formula for the dispersion)
However, it is possible to mix both ideas, i.e. specify a GLM(M) with a GLS-type correlation term (as a random effect) on the linear predictor, e.g. to account for spatial autocorrelation in GLMs, as, e.g., done in glmmPQL. 
