We can actually compute the complete asymptotic expansion of the remainder term. Introduce
$$S(x) = \sum_{k\ge 1} \left(\frac{1}{k^2}-\frac{1}{(x+k)^2}\right)$$
so that our answer is given by $$\frac{\pi^2}{6}-S(N).$$
Re-write $S(x)$ as follows:
$$S(x) = \sum_{k\ge 1} \frac{1}{k^2} \left(1-\frac{1}{(x/k+1)^2}\right).$$
The sum term is harmonic and may be evaluated by inverting its Mellin transform.
Recall the harmonic sum identity
$$\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).$
In the present case we have
$$\lambda_k = \frac{1}{k^2}, \quad \mu_k = \frac{1}{k} \quad \text{and} \quad
g(x) = 1-\frac{1}{(x+1)^2}.$$
We need the Mellin transform $h^*(s)$ of $h(x)=g(x)-1$ which is
$$\int_0^\infty -\frac{1}{(x+1)^2} x^{s-1} dx
= \left[\frac{1}{1+x} x^{s-1} \right]_0^\infty
- (s-1) \int_0^\infty \frac{1}{1+x} x^{s-2} dx.$$
The fundamental strips of these three components are $\langle 0, 2\rangle$, followed by $\langle 1,2 \rangle$ and finally, $\langle 1,2\rangle$.
The contribution in the square brackets disappears in the fundamental strip $\langle 1, 2\rangle.$
Now this last Mellin transform can be evaluated using a keyhole contour with the slot on the positive real axis, which gives
$$(1-e^{2\pi i (s-1)})\times
\int_0^\infty \frac{1}{1+x} x^{s-1} dx
= 2\pi i \times \mathrm{Res}\left(\frac{1}{1+x} x^{s-1}; x= -1\right)$$
which yields in turn
$$\int_0^\infty \frac{1}{1+x} x^{s-1} dx = 2\pi i
\frac{ e^{\pi i (s-1)}}{1-e^{2\pi i s}}
= - 2\pi i \frac{1}{e^{-\pi i s}-e^{\pi i s}}
= \frac{\pi}{\sin(\pi s)}.$$
It follows that the Mellin transform of $g(x)$ is
$$g^*(s) = - (s-1) \frac{\pi}{\sin(\pi (s-1))} = (s-1)\frac{\pi}{\sin(\pi s)}$$
with fundamental strip $\langle -1, 0\rangle$ as can be seen by expanding $g(x)$ about zero and about infinity and testing convergence of the integral.
Therefore the transform $Q(s)$ of $S(x)$ is
$$ Q(s) = (s-1)\frac{\pi}{\sin(\pi s)} \zeta(2-s)
\quad\text{because}\quad \sum_{k\ge 1}
\frac{\lambda_k}{\mu_k^s} = \sum_{k\ge 1}\frac{1}{k^2} k^s = \zeta(2-s).$$
The half plane of convergence of the zeta term is $\Re(s)<1.$
We thus obtain the Mellin inversion integral
$$ S(x) = \frac{1}{2\pi i} \int_{-1/2-i\infty}^{-1/2+i\infty} Q(s)/x^s ds$$
which we evaluate by shifting it to the right for an expansion about infinity.
We get $$\mathrm{Res}(Q(s)/x^s; s = 0) = -\frac{\pi^2}{6}$$
and
$$\mathrm{Res}(Q(s)/x^s; s = 1) = \frac{1}{x}$$
and
$$\mathrm{Res}(Q(s)/x^s; s = 2) = -\frac{1}{2x^2}.$$
For the remaining poles at $s=q$ where $q>2$ we obtain
$$\mathrm{Res}(Q(s)/x^s; s = q) = (q-1) (-1)^q \zeta(2-q)\frac{1}{x^q}
\\= (-1)^q (q-1) \zeta(-(q-2))\frac{1}{x^q}
= (-1)^{q-1} (q-1)\frac{B_{q-1}}{q-1} \frac{1}{x^q}
= (-1)^{q-1} B_{q-1}\frac{1}{x^q} .$$
Now this only contributes when $q$ is odd
so that we may simplify it to $B_{q-1}\frac{1}{x^q} .$
(Note that those zero values from the Bernoulli numbers at odd indices cancel the poles from the sine term, so it is in fact correct to let the summation for the expansion range over all residues at poles at $q>2.$)
This gives for the asymptotic expansion
$$S(x) \sim \frac{\pi^2}{6} - \frac{1}{x} + \frac{1}{2x^2}
- \sum_{q>2} \frac{B_{q-1}}{x^q}.$$
We finally have
$$\frac{\pi^2}{6} - S(N) \sim
\frac{1}{N} - \frac{1}{2N^2} + \sum_{q>2} \frac{B_{q-1}}{N^q}$$
as discovered by Ron Gordon in his excellent answer.