Let us define:
\begin{eqnarray}
J(\alpha)&:=& \int\limits_0^\infty erf\left( \frac{\alpha}{\sqrt{1+x}} - \frac{\sqrt{1+x}}{\beta} \right) e^{-\frac{x}{\gamma}} dx\\
&\underbrace{=}_{t=\sqrt{1+x}}& \int\limits_1^\infty erf\left( \frac{\alpha}{t} - \frac{t}{\beta} \right) e^{-\frac{t^2-1}{\gamma}} 2 t dt
\end{eqnarray}
Now we differentiate with respect to $\alpha$. We have:
\begin{eqnarray}
&&\frac{\partial}{\partial \alpha} J(\alpha)=\\
&& \int\limits_1^\infty \frac{4}{\sqrt{\pi}} e^{-(\frac{\alpha}{t} - \frac{t}{\beta})^2} e^{-\frac{t^2-1}{\gamma}} dt =\\
&& \frac{4}{\sqrt{\pi}} e^{2 \frac{\alpha}{\beta}+\frac{1}{\gamma}} \int\limits_1^\infty e^{-\left( \frac{\alpha^2}{t^2} + (\frac{1}{\beta^2}+\frac{1}{\gamma}) t^2\right)} dt \underbrace{=}_{\begin{array}{rrr} A&:=&\alpha^2\\B&:=&1/\beta^2+1/\gamma\end{array}}\\
&&
\frac{4}{\sqrt{\pi}} e^{2 \frac{\alpha}{\beta}+\frac{1}{\gamma}} \cdot
\left.\frac{\sqrt{\pi } \left(e^{-2 \sqrt{A} \sqrt{B}} \left(1-\text{erf}\left(\frac{\sqrt{A}}{t}-\sqrt{B} t\right)\right)+e^{2 \sqrt{A} \sqrt{B}} \left(\text{erf}\left(\frac{\sqrt{A}}{t}+\sqrt{B} t\right)-1\right)\right)}{4
\sqrt{B}} \right|_1^\infty =\\
%
&&= \frac{e^{\frac{1}{g}}}{\sqrt{\frac{1}{b^2}+\frac{1}{g}}} \left(e^{a \left(-\Delta_-\right)} \text{erf}\left(a-\sqrt{\frac{1}{b^2}+\frac{1}{g}}\right)-e^{a \left(\Delta_+\right)}
\text{erf}\left(a+\sqrt{\frac{1}{b^2}+\frac{1}{g}}\right)+e^{a \left(-\Delta_-\right)}+e^{a \left(\Delta_+\right)}\right)
\end{eqnarray}
where $\delta:=\sqrt{\frac{1}{b^2}+\frac{1}{g}}$ and $\Delta_\pm:=2 \delta \pm \frac{2}{b}
$.
Now all we need to do is to integrate over $\alpha$. A glimpse at the bottom line above suffices to realize that the anti-derivative with respect to $\alpha$ can be easily found using integration by parts. In addition it is also easy to see that $J(\infty) = \gamma$. Therefore the final result reads:
\begin{eqnarray}
\gamma- J(\alpha) =
%
\gamma\cdot (1-erf(\alpha-\frac{1}{\beta})) +
\frac{\gamma e^{\frac{1}{\gamma}} e^{-\alpha \Delta_-}}{2 \delta}
\left(\frac{1}{\beta} - \frac{\Delta_-}{2} e^{4 \alpha \delta} +\delta\right) +
\frac{\gamma e^{\frac{1}{\gamma}}}{4 \delta}
\left(
\Delta_+ e^{-\alpha \Delta_-} erf(a-\delta) + e^{\alpha \Delta_+} \Delta_- erf(a+\delta)
\right)
\end{eqnarray}
In[441]:= {b, g} = RandomReal[{0, 2}, 2, WorkingPrecision -> 50];
dd = Sqrt[1/b^2 + 1/g];
{Dm, Dp} = 2 {-1/b + dd, 1/b + dd};
a = Range[0, 2, 0.001];
f = Interpolation[
Transpose[{a,
NIntegrate[
Erf[#/Sqrt[1 + x] - Sqrt[1 + x]/b] Exp[-x/g], {x, 0,
Infinity}] & /@ a}]]
In[446]:= a = RandomReal[{0, 2}, WorkingPrecision -> 50];
f'[a]
E^(1/g)/Sqrt[
1/b^2 + 1/
g] (E^(a (2 /b - 2 Sqrt[1/b^2 + 1/g])) + E^(
a (2 /b + 2 Sqrt[1/b^2 + 1/g])) +
E^(a (2 /b - 2 Sqrt[1/b^2 + 1/g])) Erf[a - Sqrt[1/b^2 + 1/g]] -
E^(a (2 /b + 2 Sqrt[1/b^2 + 1/g])) Erf[a + Sqrt[1/b^2 + 1/g]])
g - f[a]
g (1 - Erf[a - 1/b]) + (g E^(1/g) E^(-a Dm))/(
2 dd) (1/b + E^(4 a Sqrt[1/b^2 + 1/g]) (-Dm/2) + dd) + (g E^(1/g))/(
4 dd) ((Dp) E^(-a Dm) Erf[a - dd] + E^(a Dp) (Dm) Erf[a + dd])
Out[447]= 0.648692
Out[448]= 0.648691540832040074541903643440002542062490436
Out[449]= 0.681556
Out[450]= 0.681555757815934967025234698639289684849178686