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 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

### Edit

As pointed out by Ben Bolker in a comment, PQL versus GHQ for Gaussian responses should not make a difference. However, as shown above, results of `mgcv::gam`, `mgcv::gamm` and `gamm4::gamm4` differ. Differences can be due to different optimizers used, but such differences would generally be small. The standard error differences for the parametric terms seem substantial, I don't know why, could be due to instability due to relatively small sample size. 

If you want to test for differences between the treatments, a `by` smooth might be more appriate (although it takes up more df than a factor smooth, but I honestly don't know how to use those for testing differences).

Assuming you want to take level `A` for treatment as the reference category, and check whether each of the other two levels differ, I would take the 'ordered factor' approach. This directly allows you to test whether the parametric and smooth terms differ between the different levels of `Treatment`:

    leaf_rand$oTreatment <- ordered(leaf_rand$Treatment)
    gam_model <- gamm4(N.P ~ Treatment + s(Year, k = 5) +
                           s(Year, by = oTreatment, k = 5),
                       random = ~(1|Plot),
                       data = leaf_rand, REML = TRUE)

    gam_model1 <- gamm(N.P ~ Treatment + s(Year, k = 5) +
                         s(Year, by = oTreatment, k = 5),
                       random = list(Plot=~1),
                       data = leaf_rand, REML = TRUE)

    gam_model2 <- gam(N.P ~ Treatment + s(Year, k = 5) +
                            s(Year, by = oTreatment, k = 5) +
                        s(Plot, bs = 're'), data = leaf_rand, 
                       method="REML")

    # Model summaries and plots
    summary(gam_model$gam)
    summary(gam_model1$gam)
    summary(gam_model2)

    plot(gam_model$gam)   
    plot(gam_model1$gam)   
    plot(gam_model2)              

The summaries indicate a difference between the smooth effects of `Year` between `Treatment` levels `A` and `C`. At least, the p-values are $< .05$. Only the `gamm4` model indicates a difference in the parametric effect of `Treatment` levels `A` and `C`. But point estimates are all in the same direction.

Inspecting the plots (and `edf` values) shows that the difference between the two differing levels can be described by a negative linear effect over time.

I do not see the use of comparing fit between the fitted `gam`, `gamm` and `gamm4` models. They are essentially equivalent models with identical degrees of freedom (hence the warnings), they just used different estimation and/or optimization approaches.