The authors give a simultaneous characterization of generalized inverse Gaussian (GIG) \(\mu_{p,a,b}\) and gamma distributions. The main result is as follows. Assume that the moments \(E(X^{-r-2})\), \(E(X^{-2})\), \(E(Y^r)\) and \(E(Y^{r+2})\) are finite for a fixed \(r\). If the regressions \[ E(V^{r+1}\mid U)=c_r \quad\text{and}\quad E(V^{r+2}\mid U)=c_{r+1} \] are constant, where \(V:=X^{-1}-(X+Y)^{-1}\) [see \textit{G. Letac} and \textit{J. Wesolowski}, Ann. Probab. 28, No. 3, 1371--1383 (2000; Zbl 1010.62010)], for some constants \(c_r\) and \(c_{r+1}\), then \(c_{r+1}>c_r>0\) and there exists \(a>0\) such that \(X\sim\mu_{-p,a,b}\) and \(Y\sim\Gamma(p,a)\) where \[ p=c_r/(c_{r+1}+c_r)-r>0 \quad\text{and}\quad b=1/(c_{r+1}+c_r)>0. \]