It may be interesting to observe that this sum is an example of a harmonic sum that is often used in textbooks to illustrate the use of Mellin transforms on these sums, where the trick is to compute the Mellin transform of the sum and thereafter use Mellin inversion to get an asymptotic expansion of the latter. Introduce
$$ f(x) = \sum_{k\ge 1} \left(\frac{1}{k} - \frac{1}{x+k}\right).$$
The reason this is very useful is the fact that
$$ f(n) = H_n,$$
which should be obvious. Here $H_n$ denotes the $n$-th harmonic number. So if we expand $f(x)$ at infinity we get the asymptotic expansion of harmonic numbers.
Now to do the harmonic analysis we rewrite the sum as follows:
$$ f(x) = \sum_{k\ge 1} \frac{x}{k(x+k)} =
\sum_{k\ge 1} \frac{1}{k} \frac{x}{x+k} =
\sum_{k\ge 1} \frac{1}{k} \frac{x/k}{x/k+1}.$$
Comparing this with the equation for the Mellin transform of a harmonic sum, which we recall is
$$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x); s\right) =
\left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s),$$
where $g^*(s)$ is the Mellin transform of $g(x),$ we see that in the present case
$$\lambda_k = \frac{1}{k},
\quad \mu_k = \frac{1}{k}
\quad \text{and} \quad
g(x) = \frac{x}{1+x}.$$
This gives
$$ \sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \zeta(1-s).$$
We use a keyhole contour for the Mellin transform of $g(x)$ and take into account the pole at $x=-1$, obtaining
$$ \left(1-e^{2\pi i (s-1)}\right) \int_0^\infty g(x) x^{s-1} dx =
2\pi i(-1)^s$$
which implies that
$$g^*(s) = 2\pi i\frac{e^{i\pi s}}{1-e^{2\pi i s}}
= \pi \frac{2i}{e^{-\pi i s} - e^{\pi i s}} =
-\frac{\pi}{\sin(\pi s)}.$$
We may now conclude that the Mellin transform $f^*(s)$ of $f(x)$ is given by
$$\mathfrak{M}(f(x); s) = - \zeta(1-s) \frac{\pi}{\sin(\pi s)}.$$
We are ready to apply Mellin inversion to $f^*(s)$, shifting the inversion integral, which is
$$ \frac{1}{2\pi i}\int_{-1/2-i\infty}^{-1/2+i\infty} f^*(s)/x^s\, ds,$$ to the right for an expansion at infinity.
We get
$$\operatorname{Res}(f^*(s)/x^s; s=0) = -\gamma -\log x$$
$$\operatorname{Res}(f^*(s)/x^s; s=1) = -\frac{1}{2x}$$
and from then on
$$\operatorname{Res}(f^*(s)/x^s; s=m) = - \zeta(1-m) \frac{(-1)^m}{x^m} =
\frac{(-1)^m B_m}{m}\frac{1}{x^m}$$
where the $B_m$ are Bernoulli numbers.
The conclusion is that
$$ f(x) \sim \gamma + \log x + \frac{1}{2x}
- \sum_{m\ge 2} \frac{(-1)^m B_m}{m}\frac{1}{x^m}$$
and in particular
$$ H_n \sim \gamma + \log n
+ \frac{1}{2n} - \frac{1}{12n^2} + \frac{1}{120n^4} - \cdots.$$