To add to the response by @MurrayEisenberg, this can be used to explain why the limit is what it is. Without going through all details we'll just investigate the series for Gamma[1+z] at z=0. This amounts to looking at Gamma[z] and first derivative at z=1. Recall gamma definition as
Gamma[z] == Integrate[x^(z-1) Exp[-x], {x,0,Infinity}]
So Gamma[1] is given as
In[68]:= Integrate[x^(z - 1)*Exp[-x] /. z -> 1, {x, 0, Infinity}]
Out[68]= 1
For the first derivative, we differentiate under the integral sign (omitting details as to why this is valid).
In[69]:= Integrate[D[x^(z - 1)*Exp[-x], z] /. z -> 1, {x, 0, Infinity}]
Out[69]= -EulerGamma
The rest involving the log and dividing by z is straightforward using the series of the logarithm centered at 1.
