A self-contained approach. For any $n\geq 1$ we have $\frac{1}{n}=\int_{0}^{+\infty}e^{-nt}\,dt$ and $\log(n)=\int_{0}^{+\infty}\frac{e^{-t}-e^{-nt}}{t}\,dt$ by Frullani's theorem, hence $a_n=H_n-\log(n)$ can be represented as
$$ \int_{0}^{+\infty}\frac{1-e^{-nt}}{e^t-1}-\frac{e^{-t}-e^{-nt}}{t}\,dt=\int_{0}^{+\infty}\left(\frac{1}{e^t-1}-\frac{1}{te^t}\right)\,dt-\int_{0}^{+\infty}\left(\frac{1}{e^t-1}-\frac{1}{t}\right)e^{-nt}\,dt $$
hence
$$ a_n-\gamma=\int_{0}^{+\infty}\left(\frac{1}{t}-\frac{1}{e^t-1}\right)e^{-nt}\,dt = \int_{0}^{+\infty}g(t) e^{-nt}. \tag{1}$$
We already know that the RHS goes to zero as $n\to+\infty$, and the positive function $g(t)$
- Is bounded by $\frac{1}{2}e^{-t/6}$ on the interval $(0,6]$
- Is bounded by $\frac{1}{6}$ on the interval $(6,+\infty)$
hence
$$ 0\leq a_n-\gamma \leq \frac{3}{6n+1} \tag{2}$$
for any $n\geq 2$. By applying integration by parts to the RHS of $(1)$ you may similarly show that
$$ H_n = \log(n) + \gamma + \frac{1}{2n}-\frac{1-o(1)}{12n^2}.\tag{3} $$
The same can be proved through creative telescoping or through the Euler-Maclaurin summation formula.