Here, it is useful to know the generating function ($|z|<1$)
$$ \frac{\exp\left(-\frac{xz}{1-z}\right)}{(1-z)^{\alpha+1}}
= \sum_{n=0}^\infty L^\alpha_n (x) z^n .$$ Thereby, we can evaluate the integral for all $n$ simultaneously
$$g=\sum_n z^n \int_0^\infty dt\, t^\alpha e^{-t} L_n^{\alpha+1}(t)
= \int_0^\infty dt\, t^\alpha e^{-t} \frac{\exp\left(-\frac{tz}{1-z}\right)}{(1-z)^{\alpha+2}}. $$
The integral can be brought onto the integral for the $\Gamma$ function by a change of variables
$$(1-z)^{\alpha+2} g= \int_0^\infty dt\, t^\alpha e^{ - \frac{t}{1-z} }
= (1-z)^{\alpha+1} \int_0^\infty dt \, t^\alpha e^{ -t}
=(1-z)^{\alpha+1} \Gamma(\alpha+1). $$
Thus, we have $g = \Gamma(\alpha+1)/(1-z)$ and
$$\int_0^\infty dt\, t^\alpha e^{-t} L_n^{\alpha+1}(t) = \Gamma(\alpha+1) \qquad \forall n.$$
Completeness of an orthogonal sequence of functions is a bit tricky on unbounded intervals, while it is relatively straightforward on bounded intervals. In the case of Laguerre and Hermite polynomials, there is a nice trick due to von Neumann that allows the reduction to bounded intervals.
There seems to be a bit of confusion about the interval in the statement of the question. Here's a correct statement:
For any real number $\alpha \gt -1$ the functions $\langle e^{-x/2} x^{\alpha/2} L_{n}^{(\alpha)}(x)\rangle_{n=0}^\infty$ obtained from the Laguerre polynomials $L_{n}^{(\alpha)}(x)$ are a complete orthogonal system in $L^2(0,\infty)$. The Hermite polynomials $H_n(x)$ yield the complete orthogonal system $\langle e^{-x^2/2} H_n(x)\rangle_{n=0}^\infty$ in $L^2(\mathbb{R})$.
This is proved in detail in the classic book Gábor Szegő, Orthogonal polynomials, Chapter 5. The entire chapter discusses the main properties oft the Laguerre polynomials $L^{(\alpha)}_n(x)$ for an arbitrary real number $\alpha \gt -1$ and proves their completeness in Section 5.7.
More precisely, Szegő shows in Theorem 5.7.1 on pages 108f that for fixed $\alpha \gt -1$ the functions $f_n(x) = e^{-x/2}x^{\alpha/2} x^n$ span a dense subspace of $L^2(0,\infty)$.
The first idea is to use a change of variables $y = e^{-x}$ in order to use the case of $L^2(0,1)$ where density of the span of $(\log1/y)^{\alpha/2} y^n$ is not too hard to prove (see Theorem 3.1.5).
Write a function in $L^2(0,\infty)$ as $e^{-x/2} x^{\alpha/2} f(x)$. Then we have that $(\log1/y)^{\alpha/2} f(\log(1/y)) \in L^2(0,1)$ can be approximated by functions of the form $(\log1/y)^{\alpha/2} p(y)$ where $p$ is a polynomial. Transforming back to $(0,\infty)$ this shows that
$$
\int_{0}^\infty e^{-x} x^\alpha (f(x) - p(e^{-x}))^2 \,dx \lt \varepsilon
$$
for a suitable polynomial $p$. This reduces the task to proving that for all natural $k$ there exists a polynomial $q$ such that
$$\tag{$\ast$}
\int_{0}^\infty e^{-x} x^\alpha (e^{-kx} - q(x))^2\,dx
$$
is as small as we wish.
To do this, von Neumann's trick is to use the generating function
of the Laguerre polynomials $L_{n}^{(\alpha)}(x)$
$$
(1-w)^{-\alpha-1} \exp\left(-\frac{xw}{1-w}\right) = \sum_{n=0}^\infty L_n^{(\alpha)}(x) w^n.
$$
Choosing $w = \frac{k}{k+1}$ we have $\exp\left(-\frac{xw}{1-w}\right) = \exp{(-kx)}$.
Thus, a natural choice for $q$ is $q_N(x) = (1-w)^{\alpha+1} \sum_{n=0}^N L_n^{(\alpha)}(x) w^n$ with large enough $N$. Plugging this into $(\ast)$ we obtain using the orthogonality relations
$$
\begin{align*}
\int_{0}^\infty e^{-x} x^\alpha (e^{-kx} - q_N(x))^2\,dx & = (1-w)^{2\alpha+2} \int_{0}^\infty e^{-x} x^\alpha \left(\sum_{n=N+1}^\infty L_{n}^{(\alpha)}(x) w^{n}\right)^2\,dx \\
&= (1-w)^{2\alpha+2} \Gamma(\alpha+1) \sum_{n=N+1}^\infty \binom{n+\alpha}{n} w^{2n}
\end{align*}$$
where term-wise integration is justified using an application of Cauchy-Schwarz. It remains to observe that the last expression tends to $0$ as $N \to \infty$.
Another reference discussing the case of $\alpha = 0$ nicely is Courant and Hilbert, Methods of mathematical physics, I, §9, sections 5 and 6. They discuss ordinary Laguerre and Hermite polynomials and their completeness.
Best Answer
A straightforward idea is to put the series for $J_\alpha$ directly into the integral (and exchange integration with summation, admissible because of absolute convergence). Recognizing an $n$-th derivative middleways, we reduce everything to the Rodrigues formula for $L_n^{(\alpha)}(x)$: \begin{align}\mathrm{RHS}&=x^{-\alpha}e^x\int_0^\infty y^n e^{-y}\sum_{k=0}^\infty\frac{(-1)^k(xy)^{k+\alpha}}{k!\Gamma(k+\alpha+1)}~dy\\&=x^{-\alpha}e^x\frac{\partial^n}{\partial x^n}\int_0^\infty e^{-y}\sum_{k=0}^\infty\frac{(-1)^k(xy)^{k+n+\alpha}}{k!\Gamma(k+n+\alpha+1)}~dy\\&=x^{-\alpha}e^x\frac{\partial^n}{\partial x^n}\sum_{k=0}^\infty\frac{(-1)^k x^{k+n+\alpha}}{k!\Gamma(k+n+\alpha+1)}\int_0^\infty y^{k+n+\alpha}e^{-y}~dy\\&=x^{-\alpha}e^x\frac{\partial^n}{\partial x^n}\left(x^{n+\alpha}\sum_{k=0}^\infty\frac{(-x)^k}{k!}\right)=\mathrm{LHS}.\end{align}