The Penalized Quasi Likelihood (PQL) method has been proposed to fit generalized linear mixed-effects models. The way it works is by doing a kind of a Laplace approximation in a quasi-likelihood formulation of the model. This approximation results in a transformation of the original outcome variable. The aim of the transformation is to make the transformed outcome behave like a normally distributed variable for which a linear mixed-effects model can be fitted. Hence, fitting of GLMMs can be done with software fitting linear mixed models.
This transformation works well when the original outcome vector is already close to being normally distributed. E.g., when you have Poisson data with high counts or binomial data with many trials. PQL does not work satisfactorily and results in bias when the original outcome variable is more discrete, e.g., for binary data or Poisson with low counts.
The gold standard approaches for approximating the log-likelihood function of GLMMs are the adaptive Gaussian quadrature and Monte Carlo. More information also available here.