Suppose we answer this question by computing the asymptotics of
$$\sum_{k=1}^n \frac{1}{k^3}.$$ This is actually a textbook example of
harmonic summation techniques.
Introduce
$$S(x) = \sum_{k\ge 1} \left(\frac{1}{k^3}-\frac{1}{(x+k)^3}\right)$$
so that our answer for the tail end of the series is given by
$$\zeta(3)-S(n).$$
Re-write $S(x)$ as follows:
$$S(x) = \sum_{k\ge 1} \frac{1}{k^3} \left(1-\frac{1}{(x/k+1)^3}\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^3}, \quad
\mu_k = \frac{1}{k} \quad \text{and} \quad
g(x) = 1-\frac{1}{(x+1)^3}.$$
Canceling an initial part of the series about zero of the function
being transformed only shifts the fundamental strip, so we may start
by calculating the Mellin transform $h^*(s)$ of $h(x)=1-g(x)$ which is
$$\int_0^\infty \frac{1}{(1+x)^3} x^{s-1} dx.$$
Now this 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)^3} x^{s-1} dx
= 2\pi i \times
\mathrm{Res}\left(\frac{1}{(1+x)^3} x^{s-1}; x= -1\right).$$
The branch of the logarithm used with this contour has the cut on the
positive real axis so that in a neighborhood of $x=-1$ using the
binomial series
$$x^{s-1} =
\sum_{q\ge 0}
\left.\left(\frac{d}{dx}\right)^q x^{s-1}\right|_{x=-1}
\frac{(x+1)^q}{q!}$$
and the residue is
$$\frac{1}{2!} \left.\left(\frac{d}{dx}\right)^2 x^{s-1}\right|_{x=-1}
\\= \frac{1}{2} (s-2) (s-1) (-1)^{s-3}
= \frac{1}{2} (s-2) (s-1) e^{i\pi(s-3)}
= \frac{1}{2} (s-2) (s-1) e^{i\pi(s-1)}$$
and we have for the transform that
$$h^*(s) = \int_0^\infty \frac{1}{(1+x)^3} x^{s-1} dx
= 2\pi i \frac{1}{2} (s-1) (s-2) \frac{ e^{\pi i (s-1)}}{1-e^{2\pi i s}}
\\= -\frac{1}{2} (s-1) (s-2) 2\pi i \frac{1}{e^{-\pi i s}-e^{\pi i s}}
= \frac{1}{2} (s-1) (s-2) \frac{\pi}{\sin(\pi s)}.$$
It follows that the Mellin transform of $g(x)$ is also given by
$$g^*(s) = -\frac{1}{2} (s-1) (s-2)\frac{\pi}{\sin(\pi s)}.$$
The Mellin integral that we just computed has fundamental strip
$\langle 0,3\rangle$. Canceling the first term of the expansion about
zero by adding one keeps the transform but shifts the fundamental
strip to $\langle -1,0\rangle.$
Therefore the transform $Q(s)$ of $S(x)$ is
$$ Q(s) = -\frac{1}{2} (s-1) (s-2)
\frac{\pi}{\sin(\pi s)} \zeta(3-s)
\quad\text{because}\quad \sum_{k\ge 1}
\frac{\lambda_k}{\mu_k^s} = \sum_{k\ge 1}\frac{1}{k^3} k^s
= \zeta(3-s).$$
The half plane of convergence of the zeta term is $\Re(s)<2.$
Taking into account the intersection of $\langle -1, 0\rangle$ with
the half plane of convergence 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{1}{2} \times -1 \times -2 \times (-1)^0 \times \zeta(3)
= -\zeta(3).$$
(Not to worry about the apparent sign error, since the residues will
be negated because we are shifting to the right.)
The factor $(s-1)$ cancels the pole at $s=1$ from the sine term.
The factor $(s-2)$ cancels the pole from the sine term but there is a
contribution from the zeta function term which gives
$$\mathrm{Res}(Q(s)/x^s; s = 2) =
-\frac{1}{2} \times 1 \times (-1)^2 \times -1 \times \frac{1}{x^2}
= \frac{1}{2} \frac{1}{x^2}.$$
The pole at $s=3$ is the last one to merit special treatment. We get
$$\mathrm{Res}(Q(s)/x^s; s = 3) =
-\frac{1}{2} \times 2 \times 1 \times (-1)^3 \times
\zeta(0) \times \frac{1}{x^3}
\\= \frac{1}{4} \times -2 \times \frac{1}{x^3}
= - \frac{1}{2} \frac{1}{x^3}.$$
For the remaining poles at $s=q$ where $q>3$ we obtain
$$\mathrm{Res}(Q(s)/x^s; s = q) =
-\frac{1}{2} (q-1) (q-2) (-1)^q \zeta(3-q)\frac{1}{x^q}
\\= -\frac{1}{2} (-1)^q (q-1) (q-2) \zeta(-(q-3))\frac{1}{x^q}
= \frac{1}{2} (-1)^q (q-1) (q-2) \frac{B_{q-2}}{q-2} \frac{1}{x^q}
= \frac{1}{2} (-1)^q (q-1) B_{q-2}\frac{1}{x^q} .$$
Now this only contributes when $q$ is even
so that we may simplify it to
$$\frac{1}{2} \sum_{q\ge 2} (2q-1) B_{2q-2} \frac{1}{x^{2q}}.$$
This gives for the asymptotic expansion
$$S(x) \sim \zeta(3) - \frac{1}{2x^2} + \frac{1}{2x^3}
- \frac{1}{2} \sum_{q\ge 2}
(2q-1) B_{2q-2} \frac{1}{x^{2q}}.$$
We finally have
$$\zeta(3) - S(n) \sim
\frac{1}{2n^2} - \frac{1}{2n^3}+
1/4\,{n}^{-4}-1/12\,{n}^{-6}+1/12\,{n}^{-8}-{\frac {3}{20}}\,{n}^{-10}
\\+{\frac {5}{12}}\,{n}^{-12}-{\frac {691}{420}}\,{n}^{-14}+{\frac {35}{
4}}\,{n}^{-16}-{\frac {3617}{60}}\,{n}^{-18}+{\frac {43867}{84}}\,{n}^
{-20}\\-{\frac {1222277}{220}}\,{n}^{-22}+{\frac {854513}{12}}\,{n}^{-24
}-\cdots$$
Consult this MSE link for a closely related computation.