For the sake of completeness, I will directly compute the raw moments from the density. First, under a shape/rate parametrization, the gamma distribution has density $$f_X(x) = \frac{\beta^\alpha x^{\alpha-1} e^{-\beta x}}{\Gamma(\alpha)}, \quad x > 0.$$ We will take for granted that for any choice of parameters $\alpha, \beta > 0$, we have $$\int_{x=0}^\infty f_X(x) \, dx = 1,$$ although this result is easily derived from the identity $$\int_{z=0}^\infty x^{z-1} e^{-z} \, dz = \Gamma(z).$$ Then it follows that for a positive integer $k$, $$\begin{align*} \mathrm{E}[X^k] &= \int_{x=0}^\infty x^k f_X(x) \, dx \\ &= \frac{1}{\Gamma(\alpha)} \int_{x=0}^\infty \beta^\alpha x^{\alpha+k-1} e^{-\beta x} \, dx \\ &= \frac{\Gamma(\alpha+k)}{\beta^k \Gamma(\alpha)} \int_{x=0}^\infty \frac{\beta^{\alpha+k} x^{\alpha+k-1} e^{-\beta x}}{\Gamma(\alpha+k)} \, dx \\ &= \frac{\Gamma(\alpha+k)}{\beta^k \Gamma(\alpha)}, \end{align*}$$ where in the penultimate step we observe that the integral equals $1$ because it is the integral of a gamma density with parameters $\alpha+k$ and $\beta$. For $k = 2$, we immediately obtain $\mathrm{E}[X^2] = \frac{\Gamma(\alpha+2)}{\beta^2 \Gamma(\alpha)} = \frac{(\alpha+1)\alpha}{\beta^2}.$ Another approach is via the moment generating function: $$\begin{align*} M_X(t) = \mathrm{E}[e^{tX}] &= \int_{x=0}^\infty \frac{\beta^\alpha x^{\alpha-1} e^{-\beta x + tx}}{\Gamma(\alpha)} \, dx \\ &= \frac{\beta^\alpha}{(\beta-t)^\alpha} \int_{x=0}^\infty \frac{(\beta-t)^\alpha x^{\alpha-1} e^{-(\beta-t)x}}{\Gamma(\alpha)} \, dx \\ &= \biggl(\frac{\beta}{\beta-t}\biggr)^{\!\alpha}, \quad t < \beta, \end{align*}$$ where the condition on $t$ is required for the integral to converge. We may rewrite this as $$M_X(t) = (1 - t/\beta)^{-\alpha},$$ and it follows that $$\mathrm{E}[X^k] = \left[ \frac{d^k M_X(t)}{dt^k} \right]_{t=0} = \left[(1-t/\beta)^{-\alpha-k}\right]_{t=0} \prod_{j=0}^{k-1} \frac{\alpha+j}{\beta} = \frac{\Gamma(\alpha+k)}{\beta^k \Gamma(\alpha)}.$$