Here is an approach using Mellin transforms to enrich the collection
of solutions. Write
$$S(N) = 1 + \frac{1}{N+1} + \sum_{k\ge 2} \frac{1}{\sum_{n=0}^N k^n}
= 1 + \frac{1}{N+1} + \sum_{k\ge 2} \frac{k-1}{k^{N+1}-1}
\\= 1 + \frac{1}{N+1} + \sum_{k\ge 2} \frac{k}{k^{N+1}-1}
- \sum_{k\ge 2} \frac{1}{k^{N+1}-1}.$$
There are two harmonic sums with the same base function here which we
now evaluate.
Introduce $$S_1(x, M) = \sum_{k\ge 2} \frac{1}{(xk)^M-1}
\quad\text{and}\quad S_2(x, M) = \sum_{k\ge 2} \frac{k}{(xk)^M-1}$$
so that we are interested in $S_2(1, N+1)-S_1(1, N+1).$
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 for $k\ge 2$ that
(do $S_1$ first as $S_2$ will follow)
$$\lambda_k = 1, \quad \mu_k = k
\quad \text{and} \quad
g(x) = \frac{1}{x^M-1}.$$
We need the Mellin transform $g^*(s)$ of $g(x)$ which is
$$\int_0^\infty \frac{1}{x^M-1} x^{s-1} dx.$$
We take $M = N+1 > 2$ since the sum from the beginning diverges when
$N=1.$
The fundamental strip of this Mellin transform is
$\langle 0,M\rangle.$ We will in fact use $\langle 0, M-1\rangle$
Now to evaluate this (seemingly divergent) transform we use a slice
contour with the bottom side of the slice aligned with the positive
real axis and the origin. The angle of the slice is $2\pi/M.$ The
radius of the slice is $R$ and we let it go to infinity so that we can
easily see the contribution from the arc that connects the two
straight sides vanishes by the ML bound because its length is
$\Theta(R)$ and the integrand is $$\Theta(1/R^M\times R^{(M-1)-1}) =
\Theta(1/R^2).$$
We will be integrating through the poles at $x=1$ and $x=\exp(2\pi
i/M),$ thereby picking up half the residues from these poles. The
integral along the horizontal side $\Gamma_1$ of the slice is our
transform $g^*(s)$. The contribution from the arc call it $\Gamma_2$
vanishes. Along the slanted side call it $\Gamma_3$ we put
$x=\exp(2\pi i/M)t$ obtain the following integral:
$$\int_{\Gamma_3} \frac{1}{x^M-1} x^{s-1} dx
\\= \exp(2\pi i/M)
\int_\infty^0 \frac{1}{(\exp(2\pi i/M) t)^M-1}
(\exp(2\pi i/M) t)^{s-1} dt
\\ = - \exp(2\pi i/M)
\int_0^\infty \frac{1}{t^M-1}
(\exp(2\pi i/M) t)^{s-1} dt
\\= - \exp(2\pi i \times s/M)
\int_0^\infty \frac{1}{t^M-1} t^{s-1} dt =
- \exp(2\pi i \times s/M) \times g^*(s).$$
Putting it all together we have for $g^*(s)$ that
$$g^*(s) = \left(1-e^{2\pi i \times s/M}\right)
= \pi i
\left(
\mathrm{Res}\left(\frac{1}{x^M-1} x^{s-1}; x=1\right)+
\mathrm{Res}\left(\frac{1}{x^M-1} x^{s-1}; x=\exp(2\pi i/M)\right)
\right).$$
These poles are simple and may be evaluated with a single derivative
which gives $1/(M x^{M-1})$ and produces
$$ \pi i
\left(\frac{1}{M} +
\frac{1}{M \exp(2\pi i/M)^{M-1}}
\exp(2\pi i \times (s-1)/M)\right).$$
This is
$$\frac{\pi i}{M}
\left(1 + \exp(2\pi i \times (s-1-(M-1))/M)\right)
\\= \frac{\pi i}{M}
\left(1 + \exp(2\pi i \times (s-M)/M)\right)
\\= \frac{\pi i}{M}
\left(1 + \exp(2\pi i \times s/M)\right).$$
Returning to $g^*(s)$ we finally obtain
$$g^*(s) = \frac{\pi i}{M}
\frac{1 + e^{2\pi i \times s/M}}{1-e^{2\pi i \times s/M}}
= \frac{\pi}{M}
\frac{i(e^{-\pi i \times s/M} + e^{\pi i \times s/M})}
{e^{-\pi i \times s/M}-e^{\pi i \times s/M}}
= - \frac{\pi}{M} \cot(\pi s/M).$$
It follows that the Mellin transform $Q_1(s)$ of the harmonic sum
$S_1(x,M)$ is given by
$$Q_1(s) = - \frac{\pi}{M} \cot(\pi s/M) (-1+\zeta(s))
\\ \text{because}\\
\sum_{k\ge 2} \frac{\lambda_k}{\mu_k^s} =
\sum_{k\ge 2} \frac{1}{k^s}
= -1+\zeta(s)$$
for $\Re(s) > 1.$
Similarly the Mellin transform $Q_2(s)$ of the harmonic sum
$S_2(x,M)$ is given by
$$Q_2(s) = - \frac{\pi}{M} \cot(\pi s/M) (-1+\zeta(s-1))
\\ \text{because}\\
\sum_{k\ge 2} \frac{\lambda_k}{\mu_k^s} =
\sum_{k\ge 2} k \times \frac{1}{k^s}
= -1+\zeta(s-1)$$
for $\Re(s) > 2.$
The Mellin inversion integral for the first one is
$$\frac{1}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} Q_1(s)/x^s ds$$
which we evaluate by shifting it to the right for an expansion about
infinity. For the second one we get
$$\frac{1}{2\pi i} \int_{5/2-i\infty}^{5/2+i\infty} Q_2(s)/x^s ds$$
Now note that the first pole to the right of the abscissa of
convergence is at $s=M$ and since $M\ge 3>5/2>3/2$ we may in fact join
these two inversion integrals and write
$$\frac{1}{2\pi i} \int_{5/2-i\infty}^{5/2+i\infty}
(Q_2(s)-Q_1(s))/x^s ds$$
which is
(the minus sign disappears because we are integrating clockwise in the
right half-plane)
$$\frac{1}{2\pi i} \int_{5/2-i\infty}^{5/2+i\infty}
\frac{\pi}{M} \cot(\pi s/M) (\zeta(s-1)-\zeta(s))/x^s ds.$$
Observe that
$$\mathrm{Res}\left(\cot(\pi s/M); s=qM\right) = \frac{M}{\pi}$$
with $q$ an integer.
Collecting the residues at the poles at $s = qM$ in the right half
plane and setting $x=1$ we obtain the convergent series for
$S_2(1, M)-S_1(1, M)$
$$\sum_{q\ge 1} \frac{\pi}{M} \frac{M}{\pi} (\zeta(qM-1)-\zeta(qM))
= \sum_{q\ge 1} (\zeta(qM-1)-\zeta(qM)).$$
Returning to $N$ and the sum we started with we can confirm the
earlier result that
$$S(N) = 1 + \frac{1}{N+1}
+ \sum_{q\ge 1} (\zeta(q(N+1)-1)-\zeta(q(N+1))).$$
