The differences are likely due to the different approaches functions `gamm4` and `gamm` use to approximate the likelihood. **`nlme`** (and thereby `gam` and `gamm`) uses PQL to approximate the integrands. **`lme4`** (and thereby `gamm4`) uses Gauss-Hermite quadrature.

(RE)ML estimation of GLMMs requires integrating the random effects out of the model likelihood. There is is no closed-form solution or ways to solve this analytically, so numerical methods must be used to approximate the integrals. 

From the package documentation of function `gamm4::gamm4`: "`gamm4` is based on `gamm` from package **`mgcv`**, but uses **`lme4`** rather than **`nlme`** as the underlying fitting engine via a trick due to Fabian Scheipl. `gamm4` is more robust numerically than `gamm`, and by avoiding PQL gives better performance for binary and low mean count data." 

Dimitris Rizopoulos gives a great explanation of PQL and the different ways to numerically approximate the integrals: https://stats.stackexchange.com/a/436711/173546