I found the following formula
$$\sum_{n=1}^\infty \frac{H_n}{n^q}= \left(1+\frac{q}{2} \right)\zeta(q+1)-\frac{1}{2}\sum_{k=1}^{q-2}\zeta(k+1)\zeta(q-k)$$
and it is cited that Euler proved the formula above , but how ?
Do there exist other proofs ?
Can we have a general formula for the alternating form
$$\sum_{n=1}^\infty (-1)^{n+1}\frac{H_n}{n^q}$$
$$ \begin{align} &\sum_{j=0}^k\zeta(k+2-j)\zeta(j+2)\\ &=\sum_{m=1}^\infty\sum_{n=1}^\infty\sum_{j=0}^k\frac1{m^{k+2-j}n^{j+2}}\tag{1}\\ &=(k+1)\zeta(k+4) +\sum_{\substack{m,n=1\\m\ne n}}^\infty\frac1{m^2n^2} \frac{\frac1{m^{k+1}}-\frac1{n^{k+1}}}{\frac1m-\frac1n}\tag{2}\\ &=(k+1)\zeta(k+4) +\sum_{\substack{m,n=1\\m\ne n}}^\infty\frac1{nm^{k+2}(n-m)}-\frac1{mn^{k+2}(n-m)}\tag{3}\\ &=(k+1)\zeta(k+4) +2\sum_{m=1}^\infty\sum_{n=m+1}^\infty\frac1{nm^{k+2}(n-m)}-\frac1{mn^{k+2}(n-m)}\tag{4}\\ &=(k+1)\zeta(k+4) +2\sum_{m=1}^\infty\sum_{n=1}^\infty\frac1{(n+m)m^{k+2}n}-\frac1{m(n+m)^{k+2}n}\tag{5}\\ &=(k+1)\zeta(k+4)\\ &+2\sum_{m=1}^\infty\sum_{n=1}^\infty\frac1{m^{k+3}n}-\frac1{(m+n)m^{k+3}}\\ &-2\sum_{m=1}^\infty\sum_{n=1}^\infty\frac1{m(n+m)^{k+3}}+\frac1{n(n+m)^{k+3}}\tag{6}\\ &=(k+1)\zeta(k+4) +2\sum_{m=1}^\infty\frac{H_m}{m^{k+3}} -4\sum_{n=1}^\infty\sum_{m=1}^\infty\frac1{n(n+m)^{k+3}}\tag{7}\\ &=(k+1)\zeta(k+4) +2\sum_{m=1}^\infty\frac{H_m}{m^{k+3}} -4\sum_{n=1}^\infty\sum_{m=n+1}^\infty\frac1{nm^{k+3}}\tag{8}\\ &=(k+1)\zeta(k+4) +2\sum_{m=1}^\infty\frac{H_m}{m^{k+3}} -4\sum_{n=1}^\infty\sum_{m=n}^\infty\frac1{nm^{k+3}}+4\zeta(k+4)\tag{9}\\ &=(k+5)\zeta(k+4) +2\sum_{m=1}^\infty\frac{H_m}{m^{k+3}} -4\sum_{m=1}^\infty\sum_{n=1}^m\frac1{nm^{k+3}}\tag{10}\\ &=(k+5)\zeta(k+4) +2\sum_{m=1}^\infty\frac{H_m}{m^{k+3}} -4\sum_{m=1}^\infty\frac{H_m}{m^{k+3}}\tag{11}\\ &=(k+5)\zeta(k+4) -2\sum_{m=1}^\infty\frac{H_m}{m^{k+3}}\tag{12} \end{align} $$ Letting $q=k+3$ and reindexing $j\mapsto j-1$ yields $$ \sum_{j=1}^{q-2}\zeta(q-j)\zeta(j+1) =(q+2)\zeta(q+1)-2\sum_{m=1}^\infty\frac{H_m}{m^q}\tag{13} $$ and finally $$ \sum_{m=1}^\infty\frac{H_m}{m^q} =\frac{q+2}{2}\zeta(q+1)-\frac12\sum_{j=1}^{q-2}\zeta(q-j)\zeta(j+1)\tag{14} $$
Explanation
$\hphantom{0}(1)$ expand $\zeta$
$\hphantom{0}(2)$ pull out the terms for $m=n$ and use the formula for finite geometric sums on the rest
$\hphantom{0}(3)$ simplify terms
$\hphantom{0}(4)$ utilize the symmetry of $\frac1{nm^{k+2}(n-m)}+\frac1{mn^{k+2}(m-n)}$
$\hphantom{0}(5)$ $n\mapsto n+m$ and change the order of summation
$\hphantom{0}(6)$ $\frac1{mn}=\frac1{m(m+n)}+\frac1{n(m+n)}$
$\hphantom{0}(7)$ $H_m=\sum_{n=1}^\infty\frac1n-\frac1{n+m}$ and use the symmetry of $\frac1{m(n+m)^{k+3}}+\frac1{n(n+m)^{k+3}}$
$\hphantom{0}(8)$ $m\mapsto m-n$
$\hphantom{0}(9)$ subtract and add the terms for $m=n$
$(10)$ combine $\zeta(k+4)$ and change the order of summation
$(11)$ $H_m=\sum_{n=1}^m\frac1n$
$(12)$ combine sums
Answering the first part of the question for $q$ odd we recall from the following MSE post the identity: $$ H_n = - \frac{1}{2\pi i} \int_{-1/2-i\infty}^{-1/2+i\infty} \zeta(1-s) \frac{\pi}{\sin(\pi s)}\frac{1}{n^s} ds.$$ The proof at the above cited post is sound and I will merely refer to it here since otherwise we would just include it verbatim.
This gives the formula for your sum: $$\sum_{n\ge 1} \frac{H_n}{n^q} = - \frac{1}{2\pi i} \int_{-1/2-i\infty}^{-1/2+i\infty} \zeta(1-s) \frac{\pi}{\sin(\pi s)} \zeta(q+s) ds.$$
Now shift this integral to the left to the line $\Re(s) = -1/2-(q-1),$ getting $$\sum_{n\ge 1} \frac{H_n}{n^q} = \rho_1 - \sum_{k=1}^{q-2} \zeta(1+k) (-1)^k \zeta(q-k) - \frac{1}{2\pi i} \int_{-1/2-(q-1)-i\infty}^{-1/2-(q-1)+i\infty} \zeta(1-s) \frac{\pi}{\sin(\pi s)} \zeta(q+s) ds$$ where $$\rho_1 = \operatorname{Res}\left( -\zeta(1-s) \frac{\pi}{\sin(\pi s)} \zeta(q+s); s=-(q-1)\right).$$
Make the substitution $t=s+(q-1)$ in the integral to get (not including the minus sign in front) $$ \frac{1}{2\pi i}\int_{-1/2-i\infty}^{-1/2+i\infty} \zeta(1-(t-(q-1))) \frac{\pi}{\sin(\pi (t-(q-1))} \zeta(q+t-(q-1)) dt.$$ For $q$ odd this simplifies to $$ \frac{1}{2\pi i}\int_{-1/2-i\infty}^{-1/2+i\infty} \zeta(q-t) \frac{\pi}{\sin(\pi t)} \zeta(t+1) dt.$$ Now make another substitution, namely $v=-t$, to get $$ \frac{1}{2\pi i}\int_{1/2+i\infty}^{1/2-i\infty} \zeta(q+v) \frac{\pi}{\sin(\pi v)} \zeta(1-v) dv =-\frac{1}{2\pi i}\int_{1/2-i\infty}^{1/2+i\infty} \zeta(q+v) \frac{\pi}{\sin(\pi v)} \zeta(1-v) dv$$ where the minus on the sine term cancels the one on the differential. Finally shift this integral to the line $\Re(v) = -1/2$ to obtain $$\rho_2 - \frac{1}{2\pi i}\int_{-1/2-i\infty}^{-1/2+i\infty} \zeta(q+v) \frac{\pi}{\sin(\pi v)} \zeta(1-v) dv = \rho_2 + \sum_{n\ge 1} \frac{H_n}{n^q}$$ where $$\rho_2 = \operatorname{Res}\left(- \zeta(1-v) \frac{\pi}{\sin(\pi v)} \zeta(q+v); v=0\right).$$ We have shown that $$\sum_{n\ge 1} \frac{H_n}{n^q} = \rho_1 - \sum_{k=1}^{q-2} \zeta(1+k) (-1)^k \zeta(q-k) - \left(\rho_2 + \sum_{n\ge 1} \frac{H_n}{n^q}\right).$$ This gives $$ \sum_{n\ge 1} \frac{H_n}{n^q} = \frac{1}{2} (\rho_1-\rho_2) - \frac{1}{2} \sum_{k=1}^{q-2} \zeta(1+k) (-1)^k \zeta(q-k).$$ To conclude introduce $$ W(s) = -\zeta(1-s) \frac{\pi}{\sin(\pi s)} \zeta(q+s).$$ This implies that $$ W(-s-(q-1)) = -\zeta(s+q) \frac{\pi}{\sin(\pi (-s-(q-1)))} \zeta(1-s) = - W(s)$$ because $q$ is odd. Now $$\rho_2 = \frac{1}{2\pi i} \int_{|s|=1/2} W(s) ds.$$ Put $s = -t -(q-1)$ and note that this does not change the counterclockwise orientation of the circle induced by the first integral to get $$ -\frac{1}{2\pi i} \int_{|-t-(q-1)|=1/2} W(-t-(q-1)) dt = \frac{1}{2\pi i} \int_{|-t-(q-1)|=1/2} W(t) dt = \rho_1$$ because $|-t-(q-1)|=|(-1)(t+(q-1))|=|t-(-(q-1))|.$ The conclusion is that $$ \sum_{n\ge 1} \frac{H_n}{n^q} = -\frac{1}{2} \sum_{k=1}^{q-2} \zeta(1+k) (-1)^k \zeta(q-k)$$ for $q$ odd.
Addendum. Sun Apr 27 23:57:35 CEST 2014 I don't quite see why I didn't simply evaluate the residues $\rho_1$ and $\rho_2$ as these are both easy. This does not affect the correctness of the argument.
Addendum. Sun Nov 9 23:33:24 CET 2014 In fact the equality of the two residues follows by inspection. In retrospect it appears I wanted to avoid working with the two double poles and keep everything within the limits of pen and paper.
Note that,
$\displaystyle \int_{0}^{1} x^{n-1} \mathrm{d}x = \dfrac{1}{n}$
Differentiating w.r.t. to $n$, $(p-1)$ times, we get,
$\displaystyle \dfrac{1}{n^{p}} = \dfrac{(-1)^{p-1}}{(p-1)!} \int_{0}^{1} x^{n-1} [\ln(x)]^{p-1} \mathrm{d}x$
$\displaystyle \implies \text{S} = \sum_{n=1}^{\infty} \dfrac{H_{n}}{n^{p}} = \dfrac{(-1)^{p-1}}{(p-1)!} \int_{0}^{1} [\ln(x)]^{p-1} \sum_{n=1}^{\infty} H_{n} x^{n-1} \mathrm{d}x $
Since $\displaystyle \sum_{n=1}^{\infty} H_{n} x^{n} = -\dfrac{\ln(1-x)}{1-x} $, we get,
$\displaystyle \text{S} = \dfrac{(-1)^{p}}{(p-1)!} \int_{0}^{1}\dfrac{[\ln(x)]^{p-1} \cdot \ln(1-x) }{x(1-x)} \mathrm{d}x $
Recall the Beta Function $\displaystyle \operatorname{B}(a,b) = \int_{0}^{1} x^{a-1} (1-x)^{b-1} \mathrm{d}x = \dfrac{\Gamma(a) \Gamma(b)}{\Gamma(a+b)}$
$\displaystyle \implies \text{S} = \dfrac{(-1)^{p}}{(p-1)!} \lim_{a \to 0^+} \lim_{b \to 0^+} \left(\dfrac{{\partial}^{p-1}}{\partial a^{p-1}} \left( \dfrac{\partial}{\partial b} \operatorname{B}(a,b) \right)\right) $
$\therefore \displaystyle \sum_{n=1}^{\infty} \dfrac{H_{n}}{n^{p}} = \left(1+\dfrac{p}{2} \right)\zeta(p+1)-\dfrac{1}{2}\sum_{k=1}^{p-2}\zeta(k+1)\zeta(p-k)$.
This is valid for any integer $p \geq 2$.
When $q$ is odd and greater than $1$, one can show $$ \sum_{n=1}^{\infty} \frac{H_{n}}{n^{q}} = \frac{1}{2} \sum_{k=1}^{q-2} (-1)^{k-1} \zeta(k+1) \zeta(q-k)$$
by replacing $H_{n}$ with the integral representation
$$ H_{n} = \int_{0}^{1} \frac{1-x^{n}}{1-x} \, dx \ ,$$
switching the order of integration and summation, and then repeatedly integrating by parts.
This result is also derived in Marko Riedel's answer using a different approach.
$$ \begin{align} \sum_{n=1}^{\infty} \frac{H_{n}}{n^{q}} &= \sum_{n=1}^{\infty} \frac{1}{n^{q}} \int_{0}^{1} \frac{1-x^{n}}{1-x} \, dx \\ &= \int_{0}^{1} \frac{1}{1-x} \sum_{n=1}^{\infty} \frac{1-x^{n}}{n^{q}} \, dx \\ &= \int_{0}^{1} \frac{\zeta(q)- \text{Li}_{q}(x)}{1-x} \, dx \\ &= - \Big(\zeta(q) - \text{Li}_{q}(x) \Big) \ln(1-x) \Bigg|^{1}_{0} - \int_{0}^{1} \frac{\log(1-x) \text{Li}_{q-1}(x)}{x} \, dx \\ &= -\color{#C00000} {\int_{0}^{1} \frac{\log(1-x) \text{Li}_{q-1}(x)}{x} \, dx} \\ &= \text{Li}_{2}(x) \text{Li}_{q-1}(x) \Bigg|^{1}_{0} - \int_{0}^{1} \frac{\text{Li}_{2}(x) \text{Li}_{q-2}(x)}{x} \, dx \\ &= \zeta(2) \zeta(q-1) - \int_{0}^{1} \frac{\text{Li}_{2}(x) \text{Li}_{q-2}(x)}{x} \, dx \\ &= \zeta(2) \zeta(q-1) - \text{Li}_{3}(x) \text{Li}_{q-2}(x) \Bigg|^{1}_{0} + \int_{0}^{1} \frac{\text{Li}_{3}(x)\text{Li}_{q-3}(x) }{x} \, dx \\ &= \zeta(2) \zeta(q-1) - \zeta(3) \zeta(q-2) + \int_{0}^{1} \frac{\text{Li}_{3}(x)\text{Li}_{q-3}(x) }{x} \, dx \\&= \zeta(2) \zeta(q-1) - \zeta(3) \zeta(q-2) + \zeta(4) \zeta(q-3) - \int_{0}^{1} \frac{\text{Li}_{4}(x) \text{Li}_{4-q}(x)}{x} \, dx \\ &=\zeta(2) \zeta(q-1) - \zeta(3) \zeta(q-2) + \zeta(4) \zeta(q-3) - \ldots + \zeta(q-1) \zeta(2) - \int_{0}^{1} \frac{\text{Li}_{q-1}(x) \text{Li}_{1}(x)}{x} \, dx \\ &= \sum_{k=1}^{q-2} (-1)^{k-1} \zeta(k+1) \zeta(q-k) + \color{#C00000}{\int_{0}^{1} \frac{\log(1-x) \text{Li}_{q-1}(x)}{x} \, dx} \end{align}$$
Therefore, if $q$ is odd,
$$\sum_{n=1}^{\infty} \frac{H_{n}}{n^{q}} = - \int_{0}^{1} \frac{\log(1-x) \text{Li}_{q-1}(x)}{x} \, dx = \frac{1}{2} \sum_{k=1}^{q-2} (-1)^{k-1} \zeta(k+1) \zeta(q-k).$$
Although this problem is from April 2013 I would like to take it up and try to complete the answer turning to the question
"Can we have a general formula for the alternating form?"
$$S_a(q) = \sum_{n=1}^\infty (-1)^{n+1}\frac{H_n}{n^q}$$
By inspecting the first various expressions I have made the following guess for the alternating series for even $q = 2, 4, ...$
$$S_a(q=2,4,...) = c(q)\frac{ \zeta (q+1)}{2^{q+1}}-\sum _{k=1}^{\frac{q}{2}-1} \left(1-\frac{1}{2^{q-2 k-1}}\right) \zeta (2 k+1) \zeta (q-2 k)\tag{1}$$
Here $c(q)$ are coefficients. The first 10 entries are
$$c(2,4,..,20) = \{5,59,377,2039,10229,49139,229361,1048559,4718573,20971499\}\tag{1a}$$
This sequence is not contained in https://oeis.org and I could not find a formula up to now.
For odd $q$ Mathematica returns a seemingly simple pattern
$$S_a(q=1)= \frac{\pi ^2}{12}-\frac{\log ^2(2)}{2}\tag{2a}$$
$$S_a(q=3,5,...)= \gamma \left(1-\frac{1}{2^{q-1}}\right) \zeta (q)-\;{_aF}_b^{reg}(q)\tag{2b}$$
where $\gamma$ is the Euler gamma, and ${_ aF}_b^{reg}(q)$ is the partial derivative of the regularized hypergeometric function with the parameter sets $a$ and $b$ with repect to the last parameter in $b$ taken at the argument -1.
I still need to understand this function better before posting it here. Most probably it hides a pattern similar to that of (1).
EDIT
After having completed the entry up to this point I found that the case of odd $q$ has already been treated extensively in Calculating alternating Euler sums of odd powers in March 2017.
Using these results we can easily identify the coefficients (1a) as
$$c(q) = q \left(2^q-1\right)-1$$
Divide both sides of $$\int_0^1x^{k-1}\operatorname{Li}_n(x)\ dx\overset{IBP}{=}(-1)^{n-1}\frac{H_k}{k^n}-\sum_{i=1}^{n-1}(-1)^i\frac{\zeta(n-i+1)}{k^i}$$ by $k$ then consider the summation from $k=1$ to $\infty$ we have $$\int_0^1\frac{\operatorname{Li}_n(x)}{x}\sum_{k=1}^\infty\frac{x^k}{k}\ dx=(-1)^{n-1}\sum_{k=1}^\infty\frac{H_k}{k^{n+1}}-\sum_{i=1}^{n-1}(-1)^i\zeta(n-i+1)\sum_{k=1}^\infty\frac1{k^{i+1}}$$ $$\small{-\int_0^1\frac{\operatorname{Li}_n(x)\ln(1-x)}{x}\ dx=(-1)^{n-1}\sum_{k=1}^\infty\frac{H_k}{k^{n+1}}-\sum_{i=1}^{n-1}(-1)^i\zeta(n-i+1)\zeta(i+1)}$$ Since $$-\int_0^1\frac{\operatorname{Li}_n(x)\ln(1-x)}{x}\ dx=-\sum_{k=1}^\infty\frac1{k^n}\int_0^1 x^{k-1}\ln(1-x)\ dx=\sum_{k=1}^\infty\frac{H_k}{k^{n+1}}$$ we have $$\sum_{k=1}^\infty\frac{H_k}{k^{n+1}}[1+(-1)^n]=-\sum_{i=1}^{n-1}(-1)^i\zeta(n-i+1)\zeta(i+1)$$ Let $n=2m$ $$\sum_{k=1}^\infty\frac{H_k}{k^{2m+1}}=-\frac12\sum_{i=1}^{2m-1}(-1)^i\zeta(2m+1-i)\zeta(i+1),\quad m\in\mathbb{Z}_{\ge 1}$$
We have: \begin{eqnarray} \sum\limits_{n=1}^\infty \frac{H_n}{n^q} &=& \sum\limits_{n=1}^\infty \frac{H_n}{(n+1)^q} + \zeta(q+1) \\ &=& 1/2 \left(q \zeta(q+1) - \sum\limits_{j=1}^{q-2} \zeta(j+1) \zeta(q-j) \right)+ \zeta(q+1) \end{eqnarray} where in the last line we used the result given in the answer to question Closed form expressions for harmonic sums .
Differentiate both sides of \begin{equation*} \int_0^1\frac{x^{n-1}\ln(x)\ln(1-x)}{1-x}\mathrm{d}x=\psi^{(1)}(n)[\psi(n)+\gamma]-\frac12\psi^{(2)}(n)\tag{1} \end{equation*} $(a-2)$ times w.r.t $n$ then let $n\to 1$, \begin{gather*} \int_0^1\frac{\ln^{a-1}(x)\ln(1-x)}{1-x}\mathrm{d}x=\lim_{n\to 1}\frac{d^{a-2}}{dn^{a-2}}\left(-\frac12\psi^{(2)}(n)+\psi^{(1)}(n)[\psi(n)+\gamma]\right)\\ =-\frac12\psi^{(a)}(1)+\lim_{n\to 1}\frac{d^{a-2}}{dn^{a-2}}\psi^{(1)}(n)[\psi(n)+\gamma]\\ \left\{\text{use $\frac{d^a}{dn^{a}}(f*g)=\sum_{k=0}^a \binom{a}{k} \frac{d^{a-k}}{dn^{a-k}} f* \frac{d^k}{dn^k}g$}\right\}\\ =-\frac12\psi^{(a)}(1)+\lim_{n\to 1}\sum_{k=0}^{a-2}\binom{a-2}{k}\frac{d^{a-k-2}}{dn^{a-k-2}}\psi^{(1)}(n)\frac{d^{k}}{dn^{k}}[\psi(n)+\gamma]\\ =-\frac12\psi^{(a)}(1)+\sum_{k=0}^{a-2}\binom{a-2}{k}\psi^{(a-k-1)}(1)\lim_{n\to 1}\frac{d^{k}}{dn^{k}}[\psi(n)+\gamma]\\ \{\text{separate the first term using $\psi(1)+\gamma=0$}\}\\ =-\frac12\psi^{(a)}(1)+\sum_{k=1}^{a-2}\binom{a-2}{k}\psi^{(a-k-1)}(1)\lim_{n\to 1}\frac{d^{k}}{dn^{k}}[\psi(n)+\gamma]\\ =-\frac12\psi^{(a)}(1)+\sum_{k=1}^{a-2}\binom{a-2}{k}\psi^{(a-k-1)}(1)\psi^{(k)}(1). \end{gather*}
Write $\psi^{(a)}(1)=(-1)^{a-1}a!\zeta(a+1)$, we get
$$\int_0^1\frac{\ln^{a-1}(x)\ln(1-x)}{1-x}\mathrm{d}x=\frac{1}{2}(-1)^{a}a!\zeta(a+1)$$ $$-\frac{1}{2}(-1)^{a}(a-1)!\sum_{k=1}^{a-2}\zeta(a-k)\zeta(k+1).$$
The proof completes on writing
$$\int_0^1\frac{\ln^{a-1}(x)\ln(1-x)}{1-x}\mathrm{d}x=-\sum_{n=1}^\infty H_n\int_0^1 x^n \ln^{a-1}(x)dx$$ $$=(-1)^{a}(a-1)!\sum_{n=1}^\infty\frac{H_n}{(n+1)^a}$$ $$=(-1)^{a}(a-1)!\sum_{n=1}^\infty\frac{H_{n-1}}{n^a}$$ $$=(-1)^{a}(a-1)!\left(\sum_{n=1}^\infty\frac{H_n}{n^a}-\zeta(a+1)\right)$$
Proof of $(1)$:
Differentiate both sides of the beta function: \begin{equation*} \int_0^1 x^{a-1}(1-x)^{b-1}\mathrm{d}x=\operatorname{B}(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)} \end{equation*} with respect to $a$ and $b$, \begin{equation*} \int_0^1 x^{a-1}(1-x)^{b-1}\ln(x)\ln(1-x)\mathrm{d}x=\frac{\partial^2}{\partial a\partial b}\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}. \end{equation*} Next, take the limit on both sides letting $b\to0$, \begin{gather*} \int_0^1\frac{x^{a-1}\ln(x)\ln(1-x)}{1-x}\mathrm{d}x=\lim_{ b \to 0}\frac{\partial^2}{\partial a\partial b} \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\\ =\lim_{ b \to 0}\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\left[(\psi(a)-\psi(a+b))(\psi(b)-\psi(a+b))-\psi^{(1)}(a+b) \right]\\ \left\{\text{write $\Gamma(b)=\frac{\Gamma(1+b)}{b}$ and $\psi(b)=\psi(1+b)-\frac1b$}\right\}\\ =\lim_{ b \to 0}\frac{\Gamma(a)\Gamma(1+b)}{\Gamma(a+b)}\times\\ \lim_{ b \to 0}\frac{(\psi(a)-\psi(a+b))(\psi(1+b)-\frac1b-\psi(a+b))-\psi^{(1)}(a+b)}{b}\\ \left\{\text{multiply by $b/b$ and note that $\lim_{ b \to 0}\frac{\Gamma(a)\Gamma(1+b)}{\Gamma(a+b)}=\frac{\Gamma(a)\Gamma(1)}{\Gamma(a)}=1$}\right\}\\ =\lim_{ b \to 0}\frac{(\psi(a)-\psi(a+b))(b\psi(1+b)-1-b\psi(a+b))-b\psi^{(1)}(a+b) }{b^2}\\ \{\text{now we can apply L'Hopital's rule, since we have $0/0$}\}\\ =\lim_{ b \to 0}\frac1{2b}\left\{\psi^{(1)}(a+b)\left[b\psi(a+b)-b\psi(1+b)\right]+[\psi(a)-\psi(a+b)]\right.\\ \left.[\psi(1+b)-\psi(a+b)-\psi^{(1)}(a+b)+b\psi^{(1)}(1+b)]-b\psi^{(2)}(a+b)\right\}\\ \{\text{apply L'Hopital's rule again, since we have $0/0$}\}\\ =\lim_{ b \to 0}\frac1{2}\left\{2\psi^{(1)}(a+b)[\psi(a+b)+b\psi^{(1)}(a+b)-\psi(1+b)-b\psi^{(1)}(1+b)]\right.\\ \left.+\psi^{(2)}(a+b)[b\psi(a+b)-b\psi(1+b)-1]+[\psi(a)-\psi(a+b)]\right.\\ \left.[2\psi^{(2)}(a+b)+b\psi^{(2)}(a+b)-2\psi^{(2)}(1+b)-b\psi^{(2)}(1+b)]-b\psi^{(3)}(a+b)\right\}\\ =\frac1{2}\left\{2\psi^{(1)}(a)[\psi(a)-\psi(1)]-\psi^{(2)}(a)\right\}. \end{gather*} Substitute $\psi(1)=-\gamma$ to complete the proof.
The idea of shifting gamma and digamma function then applying L'Hopital's rule is due to my friend Khalaf Ruhemi ( he is not a user on this site ).
Addendum: Following the same approach but letting $n\to1/2$, we have
$$\int_0^1\frac{\ln^{a-1}(x)\ln(1-x^2)}{1-x^2}\mathrm{d}x=2^{-a}\int_0^1\frac{\ln^{a-1}(x)\ln(1-x)}{\sqrt{x}(1-x)}\mathrm{d}x$$
\begin{gather} =-(-1)^a a!(2^{-1-a}-1)\zeta(a+1)-(-1)^a(a-1)!(2-2^{1-a})\ln(2)\zeta(a)\\ -\frac{(-1)^{a}(a-1)!}{2^{a+1}}\sum_{k=1}^{a-2} (2^{a-k}-1)(2^{k+1}-1)\zeta(a-k)\zeta(k+1). \end{gather}
Using this integral , we have
\begin{gather*} \sum_{n=1}^\infty\frac{H_n}{(2n+1)^a}=\sum_{n=1}^\infty H_n\left(\frac{(-1)^{a-1}}{(a-1)!}\int_0^1 x^{2n}\ln^{a-1}(x)\mathrm{d}x\right)\\ =\frac{(-1)^{a-1}}{(a-1)!}\int_0^1 \ln^{a-1}(x)\sum_{n=1}^\infty H_n x^{2n}\mathrm{d}x\\ =\frac{(-1)^a}{(a-1)!}\int_0^1\frac{\ln^{a-1}(x)\ln(1-x^2)}{1-x^2}\mathrm{d}x\\ =(2^{1-a}-2)\ln(2)\zeta(a)+a(1-2^{-a-1})\zeta(a+1)\nonumber\\ -\frac12\sum_{k=1}^{a-2}(2^{k+1}-1)(2^{-k}-2^{-a})\zeta(a-k)\zeta(k+1). \end{gather*}
Use the integral form of $\zeta(a)=\frac{(-1)^{a-1}}{(a-2)!}\int_0^1\frac{\ln(1-x)\ln^{a-2}(x)}{x}\mathrm{d}x$, and then interchange integration and summation signs, we get $$(-1)^aa!\sum_{k=1}^{a-1}\zeta(a-k+1)\zeta(k+1)=a(a-1)\int_0^1\int_0^1\frac{\ln(1-x)\ln(1-y)\ln^{a-2}(xy)}{xy}\textrm{d}x\textrm{d}y$$ $$\overset{x=t/y}{=}a(a-1)\int_0^1\int_0^y\frac{\ln(1-t/y)\ln(1-y)\ln^{a-2}(t)}{ty}\textrm{d}t\textrm{d}y$$ change the order of integration $$=a(a-1)\int_0^1\frac{\ln^{a-2}(t)}{t}\left(\int_t^1\frac{\ln(1-t/y)\ln(1-y)}{y}\textrm{d}y\right)\textrm{d}t$$ integrate the inner integral by parts using $\int \frac{\ln(1-t/y)}{y}\textrm{d}y=\operatorname{Li}_2(t/y)-\operatorname{Li}_2(t)$ $$=a(a-1)\int_0^1\frac{\ln^{a-2}(t)}{t}\left(\ln(1-t)(\operatorname{Li}_2(t)-\zeta(2))+\int_t^1\frac{\operatorname{Li}_2(t/y)-\operatorname{Li}_2(t)}{1-y}\textrm{d}y\right)\textrm{d}t$$ integrate by parts using the Leibniz integral rule $$=a\int_0^1\ln^{a-1}(t)\left(\frac{\ln^2(1-t)}{t}-\int_t^1\frac{\ln(1-t)-\ln(1-t/y)}{t(1-y)}\textrm{d}y\right)\textrm{d}t$$ $$=a\int_0^1\frac{\ln^{a-1}(t)}{t}\left(\ln^2(1-t)+\operatorname{Li}_2(1-t)-\zeta(2)\right)\textrm{d}t$$ $$\overset{\text{IBP}}{=}\int_0^1\ln^a(t)\left(-\frac{\ln(t)}{1-t}+\frac{2\ln(1-t)}{1-t}\right)\textrm{d}t$$ $$=(-1)^a(a+1)!\zeta(a+2)-2\sum_{n=1}^\infty H_{n-1}\int_0^1 t^{n-1}\ln^a(t)\textrm{d}t$$ $$=(-1)^a(a+1)!\zeta(a+2)-2(-1)^aa!\sum_{n=1}^\infty\frac{H_{n-1}}{n^{a+1}}$$ $$=(-1)^aa!(a+3)\zeta(a+2)-2(-1)^aa!\sum_{n=1}^\infty\frac{H_n}{n^{a+1}}$$ $$\Longrightarrow\sum_{n=1}^\infty \frac{H_n}{n^q}= \left(1+\frac{q}{2} \right)\zeta(q+1)-\frac{1}{2}\sum_{k=1}^{q-2}\zeta(k+1)\zeta(q-k).$$
Using this method, we also find:
\begin{equation*} \sum_{n=1}^\infty\frac{\overline{H}_n}{n^q}=\eta(q+1)-\frac{q}{2}\zeta(q+1)+(2-2^{1-q})\ln(2)\zeta(q)+\frac12\sum_{k=1}^{q-2}\eta(q-k)\eta(k+1), \end{equation*} and \begin{gather*} \sum_{n=1}^\infty\frac{H_n}{(2n+1)^q}=q\left(1-\frac{1}{2^{q+1}}\right)\zeta(q+1)-\left(2-2^{1-q}\right)\ln(2)\zeta(q)\nonumber\\ -\sum_{k=1}^{q-2} \left(\frac{1}{2^{q+1}}-\frac{1}{2^k}+1\right)\zeta(q-k)\zeta(k+1). \end{gather*}
For more details, check this preprint.
