$\newcommand{\al}{\alpha}
\newcommand{\de}{\delta}
\newcommand{\De}{\Delta}
\newcommand{\ep}{\varepsilon}
\newcommand{\ga}{\gamma}
\newcommand{\Ga}{\Gamma}
\newcommand{\la}{\lambda}
\newcommand{\Si}{\Sigma}
\newcommand{\thh}{\theta}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\F}{\mathcal{F}}
\newcommand{\E}{\operatorname{\mathsf E}}
\newcommand{\PP}{\operatorname{\mathsf P}}
\newcommand{\ii}[1]{\operatorname{\mathsf I}\{#1\}}
\newcommand{\tf}{\widetilde{f}}$
By straightforward manipulations, the problem can be restated as follows: Show that
\begin{equation*}
G_a:=\int_1^\infty g_a(x)\,dx \tag{1}
\end{equation*}
is increasing in $a>0$, where $g_a(x):=a^a x^{a-1} e^{-a x} /\Ga(a)$.
The key is the following construction: consider
\begin{equation*}
r_{a,b}(x):=\frac{g_b(x)}{h_{a,b}(x)},
\end{equation*}
where
\begin{equation*}
h_{a,b}(x):=g_a(\la_{a,b}(x))\la'_{a,b}(x),\quad \la_{a,b}(x):=1+\frac{b}{a}\,(x-1)-\frac12\, \Big(\frac ba - 1\Big) \ln x,
\end{equation*}
and $\la'_{a,b}(x)$ is the derivative of $\la_{a,b}(x)$ in $x$;
here everywhere $0<a\le b$ and $x>1$.
Note that $\la'_{a,b}(x)\ge\la'_{a,b}(1)= (a + b)/(2 a)>0$, $\la_{a,b}(1)=1$, and $\la_{a,b}(x)\to\infty$ as $x\to\infty$. So, by (1),
$G_a=\int_1^\infty h_{a,b}(x)\,dx$. Also, $r_{a,a}(x)=1$, and hence
\begin{equation*}
\frac{G_b-G_a}{b-a}=\int_1^\infty h_{a,b}(x)\frac{r_{a,b}(x)-r_{a,a}(x)}{b-a}\,dx,
\end{equation*}
which yields
\begin{equation*}
G'_a=\int_1^\infty g_a(x)\rho(x)\,dx,
\end{equation*}
where $G'_a$ is the derivative of $G_a$ in $a$,
\begin{equation*}
\rho(x):=\rho_a(x):=\frac{\partial r_{a,b}(x)}{\partial b}\Big|_{b=a}
=-\frac{1}{2 a x}+\frac{(a x+a-1) \ln x}{2 a x}+\ln a-\psi(a)+\frac{1}{x}-1,
\end{equation*}
and $\psi=\Ga'/\Ga$.
It suffices to show that $\rho(x)>0$. Writing $x=1+u$ (so that $u>0$), we have
\begin{equation*}
\rho'(x)2 a (1 + u)^2=a u - (a-1) \ln(1 + u)>a u - a \ln(1 + u)>0,
\end{equation*}
so that $\rho$ is increasing.
It remains to show that $\rho(1)>0$. In view of the Gauss formula
\begin{equation*}
\psi(a)=\int_0^\infty\Big(\frac{e^{-z}}{z}-\frac{e^{-az}}{1-e^{-z}}\Big)\,dz
\end{equation*}
and the elementary formula
\begin{equation*}
\ln a=\int_0^\infty\frac{e^{-z}-e^{-az}}{z}\,dz
\end{equation*}
(see e.g. Theorem 1.6.1 and formula 1.6.1, respectively, in Andrews--Askey--Roy), we see that
\begin{equation*}
\rho(1)=-\frac1{2a}+\ln a-\psi(a)
= \int_0^\infty\Big(-\frac12-\frac1{z}+\frac1{1-e^{-z}}\Big)\,e^{-az}\,dz
= \int_0^\infty R(z)\,e^{-az}\,dz,
\end{equation*}
where
\begin{equation*}
R(z):=\frac{e^z (z-2)+z+2}{2 \left(e^z-1\right) z}
\end{equation*}
is easily shown to be $>0$ for $z>0$. This completes the proof of the conjecture.
Added: Let me explain the idea of the key construction. We needed to show that $G_b>G_a$ when $0<a\le b$, with $G_a$ as in (1). In the following picture, one can see the graphs $\{(x,g_a(x))\colon 1<x<5\}$ (red) and $\{(x,g_b(x))\colon 1<x<5\}$ (green) for $a=1,b=3$.
The functions $g_a$ and $g_b$ are not directly comparable on the interval $(1,\infty)$: none of them majorizes the other one on this entire interval. The idea is then to rearrange the mass under the graph of $g_a$ so as to make a direct comparison possible. This is done by using the transformation/change of integration variable $x\mapsto\la_{a,b}(x)$, $dx\to\la'_{a,b}(x)\,dx$. The resulting function $h_{a,b}$ (whose graph is orange in the picture) is "majorized in the limit as $b\downarrow a$" by $g_b$ on the entire interval $(1,\infty)$, and the mass under the graph of $h_{a,b}$ over the interval $(1,\infty)$ is the same as that under the graph of $g_a$.
Using transportation of mass to prove inequalities is an old and very broadly used idea. One of earlier examples in this area is using the Steiner symmetrization (where segments are moved to a central position) to prove isoperimetric-type inequalities. Optimal transport by definition produces inequalities; the recent bible in this area is "Optimal Transport" by C. Villani with 846 further references, including ones to the use of transportation techniques by Marton and Talagrand to obtain concentration-of-measure inequalities.
However, I have not seen use of the mass transportation idea in contexts similar to the one here. While the idea is, in principle, very simple, it took me many failed attempts and a great deal of tying up loose ends to finally get a working transformation $\la_{a,b}$.
