semivariogram calculation - discrepancy between gls and gstat

Dear R-sig-Geo members,

a semivariogram for spatial soil data, calculated by the gls (nlme library) function gives the following estimates:

lmm.dap2b <- gls(Dap.sa ~ 1 + ID.sitio, Pd2006mm, correlation=corGaus(form=~(easting.m+northing.m)|ID.sitio, nugget=TRUE, metric="euclidean"), na.action=na.omit, method="REML")
Variogram(lmm.dap2b)

variog      dist n.pairs
1  0.3059585  2.828427     122
2  0.4040062  4.269282     125
3  0.5744688 18.110770     123
4  0.5266091 20.000000     125
....

whereas the variogram function of the gstat library gives, for similar distances (although different number of sample pairs), very different gamma values: 

 > variogram(Dap.sa ~ 1, locations = ~ easting.m + northing.m, data=Pd2006mm, cutoff=80)
    np      dist       gamma dir.hor dir.ver   id
1  219  3.225693 0.005578128       0       0 var1
2   20  6.478671 0.004254656       0       0 var1
3    3 13.513045 0.009896324       0       0 var1
4  307 19.201322 0.009708390       0       0 var1
.....  

Can anyone explain what is happening?

Thanks for any advice,
Guido Lorenz

Can anyone explain