Closed-form of $\int_0^1 B_n(x)\psi(x+1)\,dx$

The paper

V. Adamchik, Polygamma functions of negative order, J. Comp. Appl. Math. 100, (1998), 191-199)

contains the following statement (Proposition 3, Eq. (17)):

Let $n\in\mathbb{Z}_{\ge0}$ and $\Re z>0$, then \begin{align}\int_0^z x^n\psi\left(x\right) dx=&(-1)^{n-1}\zeta'(-n)+\frac{(-1)^n}{n+1}B_{n+1}H_n+\\ & \sum_{k=0}^n(-1)^kz^{n-k}{n \choose k}\left[\zeta'(-k,z)-\frac{B_{k+1}(z)H_k}{k+1}\right],\end{align} where $B_n$ and $B_n(z)$ are Bernoulli numbers and polynomials, $H_n$ are harmonic numbers, and $\zeta(s,z)$ denotes Hurwitz zeta function.

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.

This is an answer for the odd case.

Proposition. Let $k=1,2,3,\ldots$. Then

$$ \int_0^1 B_{2k+1}(x)\: \psi (x+1) \:dx=(-1)^{k+1}\frac{(2k+1)!}{(2\pi)^{2k+1}}\pi \: \zeta(2k+1)-\sum_{j=0}^{2k}\!\frac{ {{2k+1}\choose j} B_j}{2k+1-j} \quad (*) $$

The proof is here, observing that $$\begin{align} \int_0^1 B_{2k+1}(x)\: \psi (x+1) \:dx & = \int_0^{1} \frac{B_{2k+1}\left(\left\{ 1/x \right\}\right)}{x}dx. \tag1 \end{align}$$ A proof of $(1)$ may be found here.