Let $\gamma(s,x) = \int_{0}^{x}t^{s-1}e^{-t}\,dx$ and $\Gamma(s) = \int_{0}^{\infty}t^{s-1}e^{-t}\,dx$ be the incomplete gamma function and the gamma function respectively.
Claim: $\lim_{n\to\infty} \frac{\gamma(n+1,n)}{\Gamma(n+1)} = \frac{1}{2}$.
My observation: It seems like the graph of $t^{s-1}e^{-t}$ peak at $s-1$. And the graph of $t^{s-1}e^{-t}$ is (almost) symmetry about $s-1$. 
$s-1 = 30$" /> So it suggest us to think of some translation (to 0) and normalization and see what we can get from it. But the problem is, it is not a perfect symmetric shape , but tends to become more and more symmetric as $n$ increase
As you observed, the main contribution to the integral is coming from a neighbourhood of $t=n$. Thus, we change the integration variable from $t$ to $s$ by $t=ns$, to obtain $$ \gamma (n + 1,n) = \int_0^n {t^n e^{ - t} dt} = n^{n + 1} \int_0^1 {s^n e^{ - ns} ds} $$ for $n>0$. It is convenient to replace $s$ by $e^{-t}$, leading to $$ \gamma (n + 1,n) = n^{n + 1} e^{ - n} \int_0^{ + \infty } {e^{ - n(e^{ - t} + t - 1)} e^{ - t} dt} $$ for $n>0$. By the Laplace method (cf. https://en.wikipedia.org/wiki/Laplace's_method), $$ \int_0^{ + \infty } {e^{ - n(e^{ - t} + t - 1)} e^{ - t} dt} \sim \int_0^{ + \infty } {e^{ - n\frac{{t^2 }}{2}} dt} \sim \sqrt {\frac{\pi }{{2n}}} $$ as $n\to +\infty$. Hence, $$ \gamma (n + 1,n) \sim n^{n + 1/2} e^{ - n} \sqrt {\frac{\pi }{2}} $$ as $n\to +\infty$. Combining this with Stirling's formula yields the required limit. More generally, $$ \frac{{\gamma (n + 1,n + y\sqrt {2n} )}}{{\Gamma (n + 1)}} = \frac{1}{2}\operatorname{erfc}( - y) + \mathcal{O}\!\left( {\frac{1}{{\sqrt n }}} \right) $$ for bounded real $y$ as $n\to +\infty$ (cf. http://dlmf.nist.gov/8.11.E10). See also https://dlmf.nist.gov/8.12.