<s>[INST] Prove $\int_{0}^{\infty} \psi^{(2)} (1+x) \ln (x) \, dx = \zeta'(2) + \zeta(2)$ I'm looking for alternate methods of proving the following result:

$$\int_{0}^{\infty} \psi^{(2)} (1+x) \ln (x) \, dx = \frac{\pi^2}{6} \left( \gamma + \ln (2\pi)-12 \ln A +1\right) = \zeta'(2) + \zeta(2)$$
Where $\psi$ is the polygamma function, $A$ is the Glaisher-Kinkelin constant and $\gamma$ is the Euler-Mascheroni constant.

Here is how I originally proved it:
Begin with the known result $$\int_{0}^{\infty} \frac{x \ln x}{e^{2 \pi x}-1} \, dx = \frac{1}{2} \zeta'(-1)=\frac{1}{24}-\frac{1}{2} \ln A$$
then use the following well-known property of the Laplace transform
$$\int_{0}^{\infty} f(x) g(x) \, dx = \int_{0}^{\infty} \left( \mathcal{L} f\right) (y) \left( \mathcal{L}^{-1} g \right) (y) \, dy$$
Let $f(x) = \frac{x^2}{e^{2 \pi x} - 1}$ and $g(x) = \frac{\ln x}{x}$ then:
$$\left(\mathcal{L} f\right) (y) = -\frac{1}{8 \pi^3} \psi^{(2)} \left(1+\frac{y}{2\pi}\right)$$
$$\left(\mathcal{L}^{-1} g\right)(y) = -\gamma - \ln (y)$$
$$\implies \frac{1}{24} - \frac{1}{2} \ln A = \frac{1}{8\pi^3}\int_{0}^{\infty} \gamma \, \psi^{(2)} \left(1+\frac{y}{2\pi}\right) \, dy + \frac{1}{8\pi^3} \int_{0}^{\infty} \psi^{(2)} \left(1+\frac{y}{2\pi}\right) \ln (y) \, dy$$
Proceed with the substitution $y \mapsto 2\pi x \implies dy = 2 \pi\,dx$
$$\implies \frac{1+\gamma}{24} - \frac{1}{2} \ln A = \frac{1}{4\pi^2} \int_{0}^{\infty} \psi^{(2)} \left(1+x\right) \ln (2 \pi x) \, dx$$
Now separate $\ln (2\pi x) = \ln (2 \pi) + \ln (x)$ and the result quickly follows.
I was wondering if there are any alternative methods that don't rely on the initial $\frac{1}{2} \zeta'(-1)$ result.

Motivation:
I've been experimenting with the following identity:
$$\int_{0}^{\infty} \frac{x^{2n-1}}{e^{2\pi x}-1}\, dx = \frac{(-1)^{n-1} B_{2n}}{4n}\,\text{ for } n \in \mathbb{N}$$
I noticed that if one takes the partial derivative with respect to $n$ and then proceeds to take the limit as $n \to 1$, then we arrive at our $-\frac{1}{2} \zeta'(-1)$ result. This means it is possible for one to define the 'derivative' of the Bernoulli numbers at $n=1$.
More precisely:
$$B_{2}^{\prime} = \frac{1}{4} - 2 \ln A$$
is this result already known, and if so, are any other Bernoulli derivatives definable in a similar way using the Glaisher-Kinkelin constant? (I know it is possible to do it with Glaisher-Kinkelin-like constants such as the Bendersky-Adamchik constants due to their similar connection to hyperfactorials, $\zeta'(-2)$ and $\zeta(3)$ for example). I'm not sure what meaning one can give to the 'derivative' of the Bernoulli numbers, but figure 2 of this paper seems to roughly agree with my result.
 [/INST] For another solution using the complete beta function:
Using the same expansion
$$\psi^{(2)}(x) = - 2 \sum_{n=0}^{\infty} \frac{1}{(n+x)^{3}},$$
we have
$$ \begin{align} \int_{0}^{\infty} \psi^{(2)}(x+1) \ln (x) \, \mathrm dx &= -2\int_{0}^{\infty}\ln(x) \sum_{n=0}^{\infty} \frac{1}{(n+x+1)^{3}} \, \,  dx \\ &= - 2 \sum_{n=0}^{\infty} \int_{0}^{\infty} \frac{\ln (x)}{(n+x+1)^{3}} \,  dx.\\
=& - 2 \sum_{n=0}^{\infty} \frac{1}{(n+1)^3} \int_{0}^{\infty}  \frac{\ln (x)}{(1+\frac{x}{(n+1)})^{3}} \,  dx\\
=& - 2 \sum_{n=0}^{\infty} \frac{1}{(n+1)^3} \int_{0}^{\infty}  \frac{\left(\frac{d}{dt}\Big|_{t=0+}x^t\right)}{(1+\frac{x}{(n+1)})^{3}} \,  dx\\ 
=& \frac{d}{dt}\Big|_{t=0+} - 2 \sum_{n=0}^{\infty} \frac{1}{(n+1)^3} \int_{0}^{\infty}  \frac{x^t}{(1+\frac{x}{(n+1)})^{3}}dx \\
=& \frac{d}{dt}\Big|_{t=0+} - 2 \sum_{n=0}^{\infty} \frac{(n+1)^t}{(n+1)^2} \int_{0}^{\infty}  \frac{w^t}{(1+w)^{3}}dw \quad \\
=& \frac{d}{dt}\Big|_{t=0+} - 2 \sum_{n=0}^{\infty} \frac{(n+1)^t}{(n+1)^2}B(t+1,2-t)\\
=& \frac{d}{dt}\Big|_{t=0+}  \sum_{n=0}^{\infty} \frac{(n+1)^t}{(n+1)^2}\pi (t-1)t \csc(\pi t)\\
=& \sum_{n=0}^{\infty} \frac{1}{(n+1)^2} - \sum_{n=0}^{\infty} \frac{\ln(n+1)}{(n+1)^2}\\
=& \zeta(2)+\zeta'(2)
\end{align}\\$$
</s>