A possible approach may be the following one, exploiting the inverse Laplace transform.
We have:
$$ n^2\text{arccot}(n)\text{arccoth}(n)-1 = \iint_{(0,+\infty)^2}\left(\frac{\sin s}{s}\cdot\frac{\sinh t}{t}-1\right) n^2 e^{-n(s+t)}\,ds\,dt \tag{1}$$
hence:
$$ \sum_{n\geq 2}\left(n^2\text{arccot}(n)\text{arccoth}(n)-1\right)=\iint_{(0,+\infty)^2}\frac{e^{s+t} \left(1+e^{s+t}\right) (s t-\sin(s)\sinh(t))}{\left(1-e^{s+t}\right)^3 s t}\,ds\,dt\tag{2}$$
but the last integral does not look so friendly. Another chance is given by:
$$ n^2\text{arccot}(n)\text{arccoth}(n)-1 = \iint_{(0,1)^2}\left(\frac{n^2}{n^2+x^2}\cdot\frac{n^2}{n^2-y^2}-1\right)\,dx\,dy \tag{3}$$
that leads to:
$$S=\sum_{n\geq 2}\left(n^2\text{arccot}(n)\text{arccoth}(n)-1\right)=\\=\frac{3}{2}-\iint_{(0,1)^2}\left(\frac{\pi y^3\cot(\pi y)+\pi x^3\coth(\pi x)}{2(x^2+y^2)}+\frac{1}{(1+x^2)(1-y^2)}\right)\,dx\,dy \tag{4}$$
Not that appealing, but probably it can be tackled through:
$$ \cot(\pi z) = \frac{1}{\pi}\sum_{n\in\mathbb{Z}}\frac{z}{z^2-n^2},\qquad \coth(\pi z) = \frac{1}{\pi}\sum_{n\in\mathbb{Z}}\frac{z}{z^2+n^2}\tag{5} $$
that come from the logarithmic derivative of the Weierstrass product for the sine function.
That expansions can be used to derive the Taylor series of $\pi z\cot(\pi z)$ and $\pi z\coth(\pi z)$, namely:
$$ \pi z \cot(\pi z) = 1-2\sum_{n\geq 1}z^{2n}\zeta(2n),\quad \pi z \coth(\pi z) = 1-2\sum_{n\geq 1}(-1)^n z^{2n}\zeta(2n).\tag{6}$$
Since $\frac{1}{(1+x^2)(1-y^2)}=\frac{1}{x^2+y^2}\left(\frac{1}{1-y^2}-\frac{1}{1+x^2}\right)$, we also have:
$$ S = \iint_{(0,1)^2}\frac{dx\,dy}{x^2+y^2}\left(\sum_{r\geq 1}(\zeta(2r)-1)y^{2r+2}-\sum_{r\geq 1}(-1)^r(\zeta(2r)-1)x^{2r+2}\right)\tag{7} $$
By symmetry, the contributes given by even values of $r$ vanish. That gives:
$$\begin{eqnarray*} S &=& \sum_{m\geq 1}\left(\zeta(4m-2)-1\right)\iint_{(0,1)^2}\frac{x^{4m}+y^{4m}}{x^2+y^2}\,dx\,dy\\&=&\sum_{m\geq 1}\frac{\zeta(4m-2)-1}{2m}\int_{0}^{1}\frac{1+u^{4m}}{1+u^2}\,du\\&=&\sum_{m\geq 1}\frac{\zeta(4m-2)-1}{2m}\left(\frac{\pi}{2}+\int_{0}^{1}\frac{u^{4m}-1}{u^2-1}\,du\right).\tag{8}\end{eqnarray*} $$
Thanks to Mathematica, we have:
$$ \begin{eqnarray*}\sum_{m\geq 1}\frac{\zeta(4m-2)-1}{4m}&=&\int_{0}^{1}\frac{x^2 \left(4 x+\pi \left(1-x^4\right) \cot(\pi x)-\pi \left(1-x^4\right) \coth(\pi x)\right)}{4 \left(-1+x^4\right)}\,dx\\&=&-\frac{1}{24 \pi ^2}\left(10 \pi ^3+6 \pi ^2 \log\left(\frac{\pi}{4} (\coth\pi-1)\right)+6 \pi \cdot\text{Li}_2(e^{-2\pi})+3\cdot \text{Li}_3(e^{-2\pi})-3\cdot\zeta(3)\right)\end{eqnarray*} $$
that is an expected generalization of the linked question. Now we just have to deal with the missing piece.
