§1. Applying the general theory developed in §2 we find the following expression for the sum (see §3):
$$s(k) :=\sum_{n\ge1} \frac{1}{n^3(n+k)^3} = -\frac{\pi^2}{k^4}+\frac{6}{k^5}H_{k}+\frac{3}{k^4}H_{k,2}+\frac{1}{k^3}H_{k,3}\tag{1.1} $$
where $H_{k,p} = \sum_{i=1}^{k} \frac{1}{i^p}$ is the generalized harmonic number.
$(1)$ shows that your numbers $a_k$ defined more correctly with an index $k$ by the relation
$$s(k) = a_k-\frac{\pi^2}{k^4}\tag{1.2}$$
are given explicitly by
$$a_k = \frac{6}{k^5}H_{k}+\frac{3}{k^4}H_{k,2}+\frac{1}{k^3}H_{k,3}\tag{1.3}$$
The first few numbers are
$$a_k|_{k=01}^{k=10}=\{10, \frac{21}{32}, \frac{809}{5832}, \frac{2615}{55296},\frac{ 112831}{5400000}, \frac{168791}{15552000}, \frac{
17769701}{2823576000},\\
\frac{ 22201623}{5619712000},\frac{ 30715230979}{11666192832000}, \frac{
29416735711}{16003008000000}\}\tag{1.4}$$
Notice that these numbers have been given by Claude Leibovici earlier.
§2. For the derivation of even more general formulae for sums of the type
$$s(j,p,k,q) :=\sum_{n\ge 1} \frac{1}{(n+j)^p (n+k)^q}\tag{2.1}$$
we start with the expression
$$\frac{1}{(n+j)(n+k)} = \frac{1}{k-j}(\frac{1}{n+j} -\frac{1}{n+k} )=\frac{1}{k-j}\left(\left(\frac{1}{n+j}-\frac{1}{n}\right) -\left(\frac{1}{n+k}+\frac{1}{n}\right) \right)$$
which after summing over $n$ gives
$$s(j,1,k,1) = \sum_{n\ge 1} \frac{1}{(n+j)(n+k)} \\
=\frac{1}{k-j}\sum_{n\ge 1} \left(\left(\frac{1}{n+j}-\frac{1}{n}\right) -\left(\frac{1}{n+k}+\frac{1}{n}\right) \right)=\frac{1}{k-j}(H_k - H_j)\tag{2.2}$$
Here we have used the representation of the harmonic number
$$H_k = \sum_{n\ge1}(\frac{1}{n}-\frac{1}{n+k})\tag{2.3}$$
Raising the (negative) powers $p$ and $q$ can be easily done by differentiating, viz.
$$\frac{\partial}{\partial j}s(j,p,k,q) = - p\; s(j,p+1,k,q)\tag{2.4a}$$
$$\frac{\partial}{\partial k}s(j,p,k,q) = - q\; s(j,p,k,q+1)\tag{2.4b}$$
and (from $(2.3))
$$\frac{\partial}{\partial k}H_k = \sum_{n\ge1} (\frac{1}{n+k})^2=\sum_{n\ge1} \left((\frac{1}{n+k})^2-\frac{1}{n^2}\right)+ \sum_{n\ge1}\frac{1}{n^2} = - H_{k,2}+\zeta(2)\tag{2.5}$$
The higher derivatives are given recursively by
$$\frac{\partial}{\partial k}H_{k,r} = r\left(\zeta(r+1) -H_{k,r+1} \right)\tag{2.6}$$
and explicitly by
$$\frac{\partial^r}{\partial k^r}H_k =(-1)^{r+1} r!(\zeta(r+1)-H_{k,r+1}), r\ge 1\tag{2.7}$$
which follow from the definition
$$\sum_{n\ge 1}(\frac{1}{n^r} - \frac{1}{(n+k)^r} ) = H_{k,r}\tag{2.8}$$
Hence we find that the sums $s(i,p,k,q)$ can be expressed by $\zeta$-functions and (generalized) harmonic functions.
Wait! What happens when $j=k$ (a case we have tacidly excluded)? This case is just the simpler sum $s(j,p,k\to j,q) = s(j,p+q,j,0)$.
§3. Application to the specific case of the OP.
We can write
$$s(k) = \frac{1}{4} \frac{\partial^2}{\partial j^2}\frac{\partial^2}{\partial k^2}\sum_{n\ge 1} \frac{1}{(n+j)(n+k)}|_{j \to 0}\\
=\frac{1}{4} \frac{\partial^2}{\partial j^2}\frac{\partial^2}{\partial k^2}\frac{1}{k-j}(H_k - H_j)|_{j \to 0}\tag{3.1}$$
where in the last equality we have used $(2.2)$.
Now, using $(2.7)$, this can be easily converted into our main result $(1.1)$.
§4. Explicit expression for the general sum (2.1)
The shifted sum defined by $s_x(j,p,k,q) = s(j,p+1,k,q+1)$ is given by
$$\begin{align}s_x(j,p,k,q) := \sum_{n=1}^{\infty}\frac{1}{(n+j)^{p+1} (n+k)^{q+1}}\\=
+(-1)^q \binom{p+q}{q} \frac{H_j}{(j-k)^{p+q+1}}\\
+(-1)^p \binom{p+q}{q} \frac{H_k}{(k-j)^{p+q+1}}\\
+(-1)^{q} \sum _{m=1}^p \binom{p+q-m}{q} \frac{H_{j,m+1}-\zeta (m+1)}{(j-k)^{-m+p+q+1}}\\
+(-1)^{p} \sum _{n=1}^q \binom{p+q-n}{p} \frac{H_{k,n+1}-\zeta (n+1)}{(k-j)^{-n+p+q+1}}\\
\end{align}\tag{4.1}$$
Derivation
Using
$$\begin{align}& \frac{1}{(j+n)^{p+1} (k+n)^{q+1}}\\
=&\frac{(-1)^{p+q}}{p! p!} \frac{\partial ^{(p+q)}}{\partial j^{p}\, \partial k^{q}}\left(\frac{1}{(j+n) (k+n)}\right)\end{align}\tag{4.2}$$
we have
$$\begin{align}sx(j,p,k,q)
=&\frac{(-1)^{p+q}}{p! p!} \frac{\partial ^{(p+q)}}{\partial j^{p}\, \partial k^{q}}\left(\frac{H_{j}-H_{k}}{j-k}\right)\\
= &\frac{(-1)^{p+q}}{p! p!} D_{j}^{p} D_{k}^{q}\left(\frac{H_{j}-H_{k}}{j-k}\right)\\
\end{align}\tag{4.3}$$
We can write
$$\begin{align} &D_{j}^{p} D_{k}^{q}\left(\frac{H_{j}-H_{k}}{j-k}\right)
=D_{j}^{p} D_{k}^{q}\left(\frac{H_{j}}{j-k}\right)-D_{j}^{p} D_{k}^{q}\left(\frac{H_{k}}{j-k}\right)
\end{align}\tag{4.4}
$$
Carefully carrying out the derivatives with the first term
$$\begin{align} &D_{j}^{p} D_{k}^{q}\left(\frac{H_{j}}{j-k}\right)=D_{j}^{p}\left( H_{j} D_{k}^{q}\left(\frac{1}{j-k}\right)\right)\\
=&q! D_{j}^{p}\left( H_{j} \frac{1}{(j-k)^{q+1}}\right)\\
=&q!\sum _{m=0}^p \binom{p}{m} \left( D_{j}^{m} H_{j}\right) \left( D_{j}^{p-m}\frac{1}{(j-k)^{q+1}}\right)\\
=&q!\sum _{m=0}^p \binom{p}{m} \left( D_{j}^{m} H_{j}\right) \left(\frac{(-1)^{p-m} (-m+p+q)!}{q!}\frac{1}{ (j-k)^{p+q+1-m}} \right)\\
=&q!\left( H_{j}\right) \left(\frac{(-1)^{p} (p+q)!}{q!}\frac{1}{ (j-k)^{p+q+1}} \right)\\
+&q!\sum _{m=1}^p \binom{p}{m} \left( D_{j}^{m} H_{j}\right) \left(\frac{(-1)^{p-m} (-m+p+q)!}{q!}\frac{1}{ (j-k)^{p+q+1-m}} \right)\\
=& (-1)^{p} (p+q)!\frac{H_{j}}{ (j-k)^{p+q+1}}\\
+&\sum _{m=1}^p \binom{p}{m} \left( D_{j}^{m} H_{j}\right) \left(\frac{(-1)^{p-m} (-m+p+q)!}{ (j-k)^{p+q+1-m}} \right)\\
=& (-1)^{p} (p+q)!\frac{H_{j}}{ (j-k)^{p+q+1}}\\
+&\sum _{m=1}^p \binom{p}{m}(-1)^{m+1} m!\left(\zeta(m+1)-H_{j,m+1}\right) \left(\frac{(-1)^{p-m} (p+q-m)!}{ (j-k)^{p+q+1-m}} \right)
\end{align}
$$
The second term is transformed in a similar manner. Putting things togehther and simplifying gives $(4.1)$ as requested.
§5. Discussion
The general sum is composed of a transcendental part and a rational part.
Notice that this structure might be conceiled if polygamma functions are used.
The transcendental part is a linear combination of $\zeta$-functions with rational coefficients, the rational part is a similar linear combination of (generalized) harmonic numbers.
This structure is exhibited already in $(1.1)$.
The transcendental part TP of the sum
$$s_x(j,p,k,q)=\sum_{n=1}^{\infty}\frac{1}{(n+j)^{p+1} (n+k)^{q+1}}$$
for some values of $p$ and $q$ in the format $\{\{p,q\},TP(s_x)\}$ are (here $d=j-k$)
For $q=p$
$$
\begin{array}{c}
\{\{0,0\},0\} \\
\left\{\{1,1\},\frac{2 \zeta(2)}{d^2}\right\} \\
\left\{\{2,2\},-\frac{6 \zeta(2)}{d^4}\right\} \\
\left\{\{3,3\},\frac{20 \zeta(2)}{d^6}+\frac{2 \zeta(4)}{d^4}\right\} \\
\end{array}
$$
As mentioned before, for $p=q$ only even $\zeta$-functions appear.
For $q=p+1$
$$
\begin{array}{c}
\left\{\{0,1\},\frac{\zeta(2)}{d}\right\} \\
\left\{\{1,2\},\frac{\zeta(3)}{d^2}-\frac{3 \zeta(2)}{d^3}\right\} \\
\left\{\{2,3\},\frac{10 \zeta(2)}{d^5}+\frac{\zeta(4)}{d^3}-\frac{2 \zeta(3)}{d^4}\right\} \\
\left\{\{3,4\},-\frac{35 \zeta(2)}{d^7}+\frac{5 \zeta(3)}{d^6}+\frac{\zeta(5)}{d^4}-\frac{5 \zeta(4)}{d^5}\right\} \\
\end{array}
$$
For $q=p+2$
$$\begin{array}{c}
\left\{\{0,2\},\frac{\zeta(3)}{d}-\frac{\zeta(2)}{d^2}\right\} \\
\left\{\{1,3\},\frac{4 \zeta(2)}{d^4}-\frac{2 \zeta(3)}{d^3}+\frac{\zeta(4)}{d^2}\right\} \\
\left\{\{2,4\},-\frac{15 \zeta(2)}{d^6}+\frac{5 \zeta(3)}{d^5}-\frac{3 \zeta(4)}{d^4}+\frac{\zeta(5)}{d^3}\right\} \\
\left\{\{3,5\},\frac{56 \zeta(2)}{d^8}-\frac{14 \zeta(3)}{d^7}+\frac{11 \zeta(4)}{d^6}-\frac{4 \zeta(5)}{d^5}+\frac{\zeta(6)}{d^4}\right\} \\
\end{array}$$
Best Answer
This problem can be recast in terms of the famous problem of the number of ways to represent a positive integer as a sum of squares. With this perspective, we can see that the following more general statement is true for any $p > 1$ (so that each of the infinite series actually converges): $$\sum_{n=1}^{\infty} \frac{1}{(n^2)^p} + \sum_{m,n = 1}^{\infty} \frac{1}{(m^2+n^2)^p} = \left(\sum_{n=1}^{\infty} \frac{1}{n^p}\right) \left(\sum_{n=0}^{\infty} \frac{(-1)^n}{(2n+1)^p}\right).$$
The left-hand side is $$\sum_{s=1}^{\infty} \frac{n_2(s)}{s^p},$$ where $n_2(s)$ is the number of ways of representing $s$ as the sum of one or of two squares of positive integers, in which order is distinguished (i.e., $1^2 + 2^2$ is counted separately from $2^2+1^2$).
It is known that $n_2(s) = d_1(s) - d_3(s)$ (see eq. 24 in the site linked above), where $d_k(s)$ is the number of divisors of $s$ congruent to $k \bmod 4$.
The first sum on the right side of the equation has (the $p$th powers of) all positive integers as denominators and the second sum on the right has (the $p$th powers of) all odd numbers as denominators. After multiplying those sums together, then, $1/s^p$ (ignoring signs) appears on the right as many times as there are odd divisors of $s$. Each odd divisor of $s$ congruent to $1 \bmod 4$ contributes a $+1/s^p$, and each odd divisor of $s$ congruent to congruent to $3 \bmod 4$ contributes a $-1/s^p$. Thus the coefficient of $1/s^p$ on the right side is exactly $d_1(s) - d_3(s)$. Therefore, the right-hand side is also $$\sum_{s=1}^{\infty} \frac{n_2(s)}{s^p}.$$