Elementary proof:

For $s>1$, write $$\zeta(s) - \frac{1}{s-1} = \sum_{n=1}^{\infty} \left( \frac{1}{n^s} - \int_{n}^{n+1} \frac{dx}{x^s} \right).$$ I will show below that we can interchange limits and summation. Assuming this, we have $$\lim_{s \to 1^{+}} \left( \zeta(s) - \frac{1}{s-1} \right) = \sum_{n=1}^{\infty} \left( \frac{1}{n} - \int_{n}^{n+1} \frac{dx}{x} \right).$$ (It is obvious that we can interchange limits with integration by just explicitly computing the integral.)

The sum of the first $N$ terms of the above is $$\sum_{n=1}^{N} \frac{1}{n} - \int_1^{N+1} \frac{dx}{x} = \sum_{n=1}^{N} \frac{1}{n} - \log(N+1).$$ So $$\lim_{s \to 1^{+}} \left( \zeta(s) - \frac{1}{s-1} \right) = \lim_{N \to \infty} \left( \sum_{n=1}^{N} \frac{1}{n} - \log(N+1) \right)$$ and the right hand side is, by definition $\gamma$.


We now justify interchanging limits. Set $$a_n(s) = \frac{1}{n^s} - \int_{n}^{n+1} \frac{dx}{x^s}=\int_n^{n+1} \left( \frac{1}{n^s} - \frac{1}{x^s} \right) dx.$$ We are going to use the Weierstrass $M$-test, which means that we must find a convergent series $\sum M_n$ such that, for all $s$, we have $|a_n(s)| < M_n$. Actually, it is enough to do for $s$ in the interval $[1,2]$, so I will assume this from now on. I'm sure there are many ways to come up with such a boundl this is the one I found first.

We have $$a_n(s) < \int_n^{n+1} \left( \frac{1}{n^s} - \frac{1}{(n+1)^s} \right) dx = \frac{1}{n^s} - \frac{1}{(n+1)^s}=\frac{1}{n^s} \left( 1- \frac{1}{(1+1/n)^s} \right)$$ $$ < \frac{1}{n} \left( 1- \frac{1}{(1+1/n)^2} \right) < \frac{1}{n} \cdot \frac{2}{n}.$$ At the line break, we used that $1/n^s$ is maximized at the left endpoint of $[1,2]$, and $1-1/(1+1/n)^s$ is maximized at the right endpoint.

So we take $M_n = 2/n^2$ and we are done.


Assume $s > 1$. Then $$\zeta(s) = \sum_{n=1}^\infty n^{-s} =\sum_{n=1}^\infty \left( \frac{1}{\Gamma(s)} \int_0^\infty t^{s-1} \mathrm{e}^{-n t} \mathrm{d} t \right) = \frac{1}{\Gamma(s)} \int_0^\infty t^{s-1} \frac{\mathrm{e}^{-t}}{1- \mathrm{e}^{-t}} \mathrm{d} t $$ Integrability at the lower bound of integration requires $s>1$. Since $\frac{1}{1-\mathrm{e}^{-t}} = t^{-1} + \frac{1}{2} + \mathcal{o}(t)$, we subtract the main term and split it off. We can do this for $s>1$ as both integrals converge: $$ \begin{eqnarray} \zeta(s) &=& \frac{1}{\Gamma(s)} \int_0^\infty t^{s-1} \mathrm{e}^{-t} \left( \frac{1}{1-\mathrm{e}^{-t}} - \frac{1}{t} \right) \mathrm{d} t + \frac{1}{\Gamma(s)} \int_0^\infty t^{s-2} \mathrm{e}^{-t} \mathrm{d} t \\ &=& \frac{1}{\Gamma(s)} \int_0^\infty t^{s-1} \mathrm{e}^{-t} \left( \frac{1}{1-\mathrm{e}^{-t}} - \frac{1}{t} \right) \mathrm{d} t + \frac{\Gamma(s-1)}{\Gamma(s)} \\ &=& \frac{1}{\Gamma(s)} \int_0^\infty t^{s-1} \mathrm{e}^{-t} \left( \frac{1}{1-\mathrm{e}^{-t}} - \frac{1}{t} \right) \mathrm{d} t + \frac{1}{s-1} \end{eqnarray} $$ The remaining integral is finite at $s=1$, and gives the free term in the expansion you seek. Accidentally, it also gives the integral representation for the Euler-Mascheroni constant: $$ \gamma = \int_0^\infty \mathrm{e}^{-t} \left( \frac{1}{1-\mathrm{e}^{-t}} - \frac{1}{t} \right) \mathrm{d} t $$


Ok, first note that $\displaystyle\gamma = \lim (\,\sum_1^n \frac{1}{k} - \log k) = \lim(\,\sum_1^n \frac{1}{k} - \log (k + 1))$.

For $x>1$, $\displaystyle\frac{1}{x-1} = \int_1^\infty \frac{\mathrm{d} t}{t^x} = \sum_{n = 1}^\infty \int_n^{n = 1} \frac{\mathrm{d} t}{t^x}$

In addition, $\displaystyle \log (n + 1) = \int_1^{n+1} \frac{\mathrm{d}t}{t} = \sum_{k = 1}^n \int_k^{k+1} \frac{\mathrm{d}t}{t}$

So $\displaystyle\lim_{x \to 1^+} \sum_1^\infty \left( \frac{1}{n^x} - \frac{1}{x^n}\right) = \lim_{x \to 1^+} \sum_1^\infty \left(\frac{1}{n^x} - \int_n^{n+1} \frac{\mathrm{d}t}{t^x}\right)$ on the one hand and

$\displaystyle \gamma = \lim \left( \sum \frac{1}{n} - \log n \right) = \sum \left( \frac{1}{n} - \int_n^{n+1} \frac{\mathrm{d}t}{t} \right)$ on the other hand.

The second is a termwise limit of the former, and so they are equal. An application of the Weierstrass M-test gives uniform convergence of $\displaystyle\lim_{x \to 1^+} \sum_1^\infty \left( \frac{1}{n^x} - \frac{1}{x^n}\right)$, so that we can rearrange the order and we get

$$\lim_{x \to 1^+} \left( \zeta(x) - \frac{1}{x-1} \right) = \gamma$$