One fruitful way to approximate the root(s) of the Taylor polynomials is to first approximate the polynomials themselves then approximate the roots of the approximation. If you're careful you can get good error bounds between the approximate zeros and the actual zeros.
A simple way to get an integral representation for the Taylor polynomials is to use the integral form for the Taylor remainder:
$$
e^x = s_n(x) + \frac{1}{n!} \int_0^x e^t (x-t)^n\,dt,
\tag{1}
$$
where
$$
s_n(x) := \sum_{k=0}^{n} \frac{x^k}{k!}.
$$
Making the substitution $t = xs$ in the integral then yields
$$
e^x = s_n(x) + \frac{x^{n+1}}{n!} \int_0^1 e^{xs}(1-s)^n\,ds.
\tag{2}
$$
Inside the integral, if we were to send $n \to \infty$, we might notice that interesting things could happen when $x$ is proportional to $n$. Then the factors $e^{xs}$ and $(1-s)^n$ will both have an "$n$" in the exponent (and hence balance in a sense), and we can expect that the integral will tend to different values depending on the constant of proportionality. This suggests that we replace $x$ by $nx$ in the equation and thus consider
$$
\begin{align}
e^{nx} &= s_n(nx) + \frac{(nx)^{n+1}}{n!} \int_0^1 e^{nxs}(1-s)^n\,ds \\
&= s_n(nx) + \frac{(nx)^{n+1}}{n!} \int_0^1 [e^{xs}(1-s)]^n\,ds \\
&= s_n(nx) + \frac{(nx)^{n+1}}{n!} \int_0^1 \exp\{n[xs + \log(1-s)]\}\,ds.
\tag{3}
\end{align}
$$
For $x < 1$ the exponent
$$
\varphi(s) = xs + \log(1-s)
$$
is strictly decreasing on the interval $(0,1)$ and near its maximum at $s=0$ we have
$$
\varphi(s) = (x-1)s + O\!\left(s^2\right).
$$
Following the usual steps of the Laplace method we can therefore deduce that
$$
\int_0^1 \exp\{n[xs + \log(1-s)]\}\,ds = \int_0^\infty e^{n[(x-1)s]}\,ds + O\!\left(\frac{1}{n^2}\right),
\tag{4}
$$
where the constant in the error term is uniform as long as $x$ remains bounded away from $1$. Of course we can calculate the remaining integral exactly:
$$
\int_0^\infty e^{n[(x-1)s]}\,ds = \frac{1}{n} \frac{1}{1-x},
$$
and substituting these in $(3)$ yields
$$
e^{nx} = s_n(nx) + \frac{n^n}{n!} \frac{x^{n+1}}{1-x} \left[1 + O\!\left(\frac{1}{n}\right)\right].
\tag{5}
$$
A little help from Stirling's formula
$$
\frac{n^n}{n!} = \frac{e^n}{\sqrt{2\pi n}} \left[ 1 + O\!\left(\frac{1}{n}\right)\right]
$$
allows us to simplify this a bit and we end up with
$$
e^{nx} = s_n(nx) + \frac{e^n x^{n+1}}{\sqrt{2\pi n}(1-x)} \left[1 + O\!\left(\frac{1}{n}\right)\right].
\tag{6}
$$
Finally, if $x$ is a zero of $s_n(nx)$ then we may conclude that it satisfies
$$
e^{nx} = \frac{e^n x^{n+1}}{\sqrt{2\pi n}(1-x)} \left[1 + O\!\left(\frac{1}{n}\right)\right]
$$
or, after a bit of rearranging,
$$
\left(x e^{1-x}\right)^{-n} = \frac{x}{\sqrt{2\pi n}(1-x)} \left[1 + O\!\left(\frac{1}{n}\right)\right].
\tag{7}
$$
This is essentially equation $(1)$ in my answer here and can be used as a starting point to approximate the real zero $x$ (which tends to $-W(1/e)$ as $n \to \infty$ with $n$ odd) using the methods indicated there.
If I recall correctly you should be able to use $(7)$ to calculate the terms in the asymptotic series for the root up to order $O(n^{-2})$ (which is the appoximation I gave in my linked answer). To calculate more terms in the series you can be more careful with the application of the Laplace method in equation $(4)$ above; if you calculate the full asymptotic expansion for the integral in $(4)$ you should, in principle, be able to calculate any term in the asymptotic expansion of the root $x$.
I am not aware of a closed form for the coefficients in the asymptotic expansion of the root though I suppose one could be obtained using something like the Lagrange inversion formula.
