Derivation
We have to calculate the integral
$$f_{2}(a) =\int_0^{\infty } e^{-a q} \Gamma (0,q)^2 \, dq \tag{1}$$
with $a\ge0$. Here
$$\Gamma (0,q) = \int_{q}^{\infty} \frac{e^{-t}}{t}\,dt$$
is the incomplete Gamma function.
We start with partial integration writing
$$a e^{-a q} \Gamma (0,q)^2= -\frac{\partial \left(e^{-a q} \Gamma (0,q)^2\right)}{\partial q}-\frac{2 e^{-(a+1) q} \Gamma (0,q)}{q}\tag{2}$$
Here we have used
$$\frac{\partial \Gamma (0,q)^2}{\partial q}=-\frac{2 e^{-q} \Gamma (0,q)}{q}\tag{3}$$
Before we integrate $(2)$, in order to avoid the divergence at $q=0$, we differentiate $(2)$ with respect to $a$, leading to
$$e^{-a q} \Gamma (0,q)^2-a q e^{-a q} \Gamma (0,q)^2=2 e^{-(a+1) q} \Gamma (0,q)+\frac{\partial \left(q e^{-a q} \Gamma (0,q)^2\right)}{\partial q}\tag{4}$$
Now we integrate $(4)$ over $q$, which gives, identifying $f_{2}(a)$ and its derivative
$$f_{2}(a) +a \frac{\partial f_{2}(a)}{\partial a}=\frac{\log (a+2)}{a+1}\tag{5}$$
Here we have used that one integral can be done explicitly
$$\int_0^{\infty } 2 e^{-(a+1) q} \Gamma (0,q) \, dq=\frac{\log (a+2)}{a+1}\tag{6}$$
and that the integrated term
$$i = q e^{-a q} \Gamma (0,q)^2\tag{7}$$
vanishes at the boundaries $q\to0$ and $q\to\infty$. Notice that since $\Gamma (0,q)\sim \log(q)$ the factor $q$ is important for $i$ to vanish in the limit $q\to0$.
Hence with $(5)$ we have derived a ODE for the function to be found.
It is a linear inhomogeneous ODE of first order.
The solution is found (with Mathematica) to be
$$f_{2}(a) = \frac{c}{a}-\frac{2 \operatorname{Li}_2(-a-1)}{a}$$
The constant of integration $c$ can be determined by the condition that $f_{2}(a) \lt \infty$ for $a\to 0$.
A sufficient condition for this is that
$$c = 2 \operatorname{Li}_2(-1)= -\frac{\pi^2}{6}$$
This gives
$$f_{2}(a) =-\frac{1}{a}\left(2 \operatorname{Li}_2(-a-1)+\frac{\pi ^2}{6}\right)$$
Taking the limit $a\to 0$ gives
$$\lim_{a\to 0} \, f_{2}(a) = \log(4)$$
and this result can also be drived from $(1)$ with a simplified version of this derivation. This completes the derivation.
Extension
With the same method I have also calculated the integral with the third power of $\Gamma(0,q)$
$$f_{3}(a) =\int_0^{\infty } e^{-a q} \Gamma (0,q)^3 \, dq \tag{e1}$$
The ODE here is
$$f_{3}(a)+a \frac{\partial f_{3}(a)}{\partial a}=3 f_{2}(a+1)\tag{e2}$$
The initial value at $a=0$ is
$$f_{3}(a\to0) = \int_0^{\infty } \Gamma (0,q)^3 \, dq = \left(-6 \operatorname{Li}_2(-2)-\frac{\pi ^2}{2}\right)\tag{e3}$$
The lengthy result which I could not simplify further but which I checked numerically is
$$f_{3}(a) =\frac{-1}{a}\left(6 \operatorname{Li}_3(-a-2)-6 \operatorname{Li}_3(-a-1)+6 \operatorname{Li}_3\left(-\frac{2}{a+1}\right)+6 \operatorname{Li}_3\left(-\frac{a+3}{a+1}\right)+6 \operatorname{Li}_3\left(\frac{2}{a+3}\right)+6 \operatorname{Li}_3\left(\frac{a+2}{a+3}\right)+6 \operatorname{Li}_2(-a-2) \log (a+1)+6 \operatorname{Li}_2(-a-1) \log (a+1)-6 \operatorname{Li}_2(-a-2) \log (a+3)+6 \operatorname{Li}_2\left(-\frac{2}{a+1}\right) \log \left(\frac{a+1}{a+3}\right)+6 \operatorname{Li}_2\left(-\frac{a+3}{a+1}\right) \log \left(\frac{a+1}{a+3}\right)+4 \log ^3(a+1)-2 \log ^3(a+3)-3 \log (2) \log ^2(a+1)-9 \log (a+3) \log ^2(a+1)+6 \log ^2(a+3) \log (a+1)-3 \log (2) \log ^2(a+3)-3 \log (a+2) \log ^2(a+3)+\frac{1}{2} \pi ^2 \log (a+1)+6 \log (2) \log (a+3) \log (a+1)+6 \log (a+2) \log (a+3) \log (a+1)-6 \operatorname{Li}_3(-3)-12 \operatorname{Li}_3(-2)-12 \operatorname{Li}_3\left(\frac{2}{3}\right)+6 \operatorname{Li}_2(-3) \log (3)+12 \operatorname{Li}_2(-2) \log (3)-\frac{9 \zeta (3)}{2}+2 \log ^2(3) (3 \log (2)+\log (3))\right)\tag{e4}$$
In the limit $a\to 0$ we get
$$\lim_{a\to 0} \, f_{3}(a)=\int_0^{\infty }\Gamma (0,q)^3 \, dq\\=-5 \operatorname{Li}_2(-2)+\operatorname{Li}_2\left(\frac{2}{3}\right)-\frac{\pi ^2}{2}\\+\frac{\log ^2(3)}{2}-8 \log (2) \log (3)+4 \log (4) \log (3)\simeq 3.68568\tag{e5}$$
Comparing this with the much simpler expression $(e3)$ leaves some hope that the expression for $f_{3}(a)$ could be simplified appreciably.
The general case
Luckily, the general integral
$$f_{n}(a) =\int_0^{\infty } e^{-a q} \Gamma (0,q)^n \, dq $$
can be found for any $n=0,1,2,...$.
In fact, writing the integrand as
$$e^{-a q} \Gamma (0,q)^n=\frac{\partial \left(q e^{-a q} \Gamma (0,q)^n\right)}{\partial q}-q \frac{\partial \left(e^{-a q} \Gamma (0,q)^n\right)}{\partial q}$$
or, using
$$\frac{\partial \left(e^{-a q} \Gamma (0,q)^n\right)}{\partial q}=-\frac{n e^{-a q-q} \Gamma (0,q)^{n-1}}{q}-a e^{-a q} \Gamma (0,q)^n$$
as
$$e^{-a q} \Gamma (0,q)^n=n e^{-(a+1) q} \Gamma (0,q)^{n-1}+a q e^{-a q} \Gamma (0,q)^n+\frac{\partial \left(q e^{-a q} \Gamma (0,q)^n\right)}{\partial q}$$
and integrating over $q$ from $q=0$ to $q=\infty$, observing that the integrated part
$$i = q e^{-a q} \Gamma (0,q)^n$$
vanishes at both ends, we find the following ODE for $f_{n}(a)$
$$f_{n}(a) + a f'_{n}(a) = n f_{n-1}(1+a)\tag{g1}$$
Defining
$$g_{n}(a) = \frac{1}{a} f_{n}(a)\tag{g2}$$
we have, even simpler,
$$g'_{n}(a) = \frac{n}{a+1}g_{n-1}(a+1)\tag{g3}$$
with the initial condition
$$g_{n}(0)=0\tag{g4a}$$
and the "boundary" condition
$$g_{0}(a) = 1\tag{g4b}$$
Summarizing, the equations $g(3)$ and $g(4)$ form a complete recursive set of equations for the solutions $g_{n}(a)$ and, using $(g2)$, for the general integral $f_{n}(a)$.
QED.
Discussion
1) @Le Blanc observed in a solution that the problem can be naturally formulated in terms of the Laplace transform. Unfortunatey, this solution was removed.
But we can now reverse the logic and say that we have calculated the Laplace transform of the powers of the incomplete Gamma function.
2) It must be pointed out that the general scheme does not tell us much about possible closed expressions. On the contrary: already for $n=4$ we meet integrals like these, where I'm not aware of closed solutions:
$$i_{1}= \int \frac{\log(1+a) \log(2+a) \log(3+a)}{a+1}\,da$$
$$i_{2}=\int \frac{\operatorname{Li}_3(a+3)}{a+1} \, da$$
$$i_{3}=\int \frac{\text{Li}_3\left(-\frac{a+3}{a+1}\right)}{a+1} \, da$$
