I have a question about higher-order asymptotics of generalized hypergeometric functions. According to https://dlmf.nist.gov/15.4 the following is well known: $$ _2F_1(a,b;a+b;z)\sim -\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\log(1-z),\ \ z\rightarrow1^{-}. $$ My collaborator was able to coax Wolfram Mathematica into giving a higher-order estimate $$ _2F_1(a,b;a+b;z)\sim -\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\left[\log(1-z)+\psi(a)+\psi(b)+2\gamma\right]+\mathcal{O}((1-z)\log(1-z)) $$ as $z\rightarrow 1^{-}$ for $a$ and $b$ real and positive and where $\psi(z)$ is the digamma function and $\gamma$ is the Euler–Mascheroni constant.
Can anybody provide me with any reference or a hint of why the above Mathematica estimate would be correct?