Is there a closed-form of the following integral?
$$I_n = \int_0^1 B_n(x)\psi(x+1)\,dx,$$
where $B_n(x)$ are the Bernoulli polynomials and $\psi(x)$ is the digamma function.
The motivation of the problem was this question. In this answer achille hui showed that for all $n>0$
$$\int_0^1 B_n\left(\left\{ \frac1x \right\}\right) \frac{dx}{x} = -\int_0^1 B_n(u) \psi(1+u) du,$$
where $\{x\}$ is the fractional part of $x$.
By looking at specific $n$ values we could observe interesting patterns.
$$\begin{align} I_0 & = 0\\ I_1 & = 1-\frac 12 \ln(2\pi) \\ I_2 & = -\frac 12 + 2 \ln A \\ I_3 & = \frac{1}{12} - \frac{3}{4\pi^2}\zeta(3) = \frac{1}{12} - \frac{1}{8}\frac{\zeta(3)}{\zeta(2)}\\ I_4 & = \frac{1}{45} - 4\zeta'(-3)\\ I_5 & = -\frac{13}{360} + \frac{15}{4\pi^4}\zeta(5) = -\frac{13}{360} + \frac{1}{24}\frac{\zeta(5)}{\zeta(4)}\\ I_6 & = -\frac{1}{252}-6\zeta'(-5)\\ I_7 & = \frac{47}{1260} - \frac{315}{8\pi^6}\zeta(7) = \frac{47}{1260} - \frac{1}{24}\frac{\zeta(7)}{\zeta(6)}\\ I_8 & = -\frac{8}{1575}-8\zeta'(-7)\\ I_9 & = -\frac{1703}{25200} -\frac{2835}{4\pi^8}\zeta(9) = -\frac{1703}{25200} -\frac{3}{40}\frac{\zeta(9)}{\zeta(8)}\\ I_{10} & = \frac{2461}{83160} - 10\zeta'(-9)\\ \dots \end{align}$$
where $A$ is the Glaisher–Kinkelin constant, $\zeta(3)$ is the Apéry's constant and in general $\zeta$ is the Riemann zeta function and $\zeta'$ is its derivative.
The paper
contains the following statement (Proposition 3, Eq. (17)):
Specializing to $z=1$, we get $$\int_0^1 x^n\psi\left(x\right) dx=\sum_{k=0}^{n-1}(-1)^k{n \choose k}\left[\zeta'(-k)-\frac{B_{k+1}H_k}{k+1}\right].\tag{1}$$ Using (1) and the recursion relation $\psi(x+1)-\psi(x)=x^{-1}$, the integral $\int_0^1 P\left(x\right)\psi\left(x+1\right) dx$ can be computed for any polynomial $P(x)$.
Edit: There remains the task of explaining why Bernoulli polynomials provide a better basis than monomials $x^n$ (instead of linear combination of $\zeta'(-k)$ we get only one of these values). I think this could be understood after careful reading of the paper.