I developed the following:
Consider $$\frac{t}{e^t-1}=\sum_{k=0}^\infty\frac{B_k}{k!}t^k,$$ where $B_k$ are the Bernoulli numbers, then $$\frac{e^t}{e^{e^t}-1}=\sum_{k=0}^\infty\frac{B_k}{k!}e^{tk}=\sum_{k=1}^\infty\frac{B_k}{k!}\sum_{n=0}^\infty\frac{(tk)^n}{n!}=\sum_{n=0}^\infty\frac{t^n}{n!}\sum_{k=1}^\infty\frac{k^n}{k!}B_k.$$ Using a CAS to compute the Taylor series of the LHS suggests the ansatz $$\sum_{k=1}^\infty\frac{k^n}{k!}B_k=\frac{P_n(e)}{(e-1)^{n+1}},$$ where $P_n(x)$ are $n$-th degree polynomials. With this direction in mind, using Cauchy's integral formula with the generating function for the Bernoulli numbers, $$\sum_{k=0}^\infty\frac{k^n}{k!}B_k=\frac{1}{2\pi i}\sum_{k=1}^\infty k^n\oint\frac{z}{e^z-1}\frac{dz}{z^{k+1}}=\frac{1}{2\pi i}\oint\frac{\text{Li}_{-n}(1/z)}{e^z-1}.$$ Applying the closed form for the polylogarithm of negative index, $$\frac{1}{2\pi i}\sum_{k=0}^{n-1}A_{n,k}\oint\frac{1}{e^z-1}\frac{z^{k+1}}{(z-1)^{n+1}},$$ where $A_{n,k}$ are the Eulerian numbers. Using a circular contour to include the pole at $z=1$ but avoiding the poles at $\pm 2\pi i$, we obtain $$\frac{1}{n!}\left.\frac{d^n}{dz^n}\frac{zA_n(z)}{e^z-1}\right|_{z=1}=\frac{1}{n!}\sum_{k=0}^{n-1}A_{n,k}\left.\frac{d^n}{dz^n}\frac{z^{k+1}}{e^z-1}\right|_{z=1}.$$ Rewriting, we can express as $$\frac{d^n}{dz^n}\frac{e^{-z}z^{k+1}}{1-e^{-z}}=\frac{d^n}{dz^n}\sum_{\ell=1}^\infty e^{-z\ell}z^{k+1}=(-1)^{n-k+1}(k+1)!\sum_{\ell=1}^\infty e^{-\ell} \ell^{n-k-1} L_{k+1}^{n-k-1}(\ell),$$ where $L_a^b(x)$ are the generalised Laguerre polynomials. Hence, $$\sum_{k=0}^\infty\frac{k^n}{k!}B_k=\frac{(-1)^{n+1}}{n!}\sum_{k=0}^{n-1}A_{n,k}(-1)^k(k+1)!\sum_{\ell=1}^\infty e^{-\ell}\ell^{n-k-1}L_{k+1}^{n-k-1}(\ell).$$ Now $$\frac{(-1)^{n+1}}{n!}\sum_{k=0}^{n-1}A_{n,k}(-1)^k(k+1)!\sum_{\ell=1}^\infty e^{-\ell}\ell^{n-k-1}\sum_{i=0}^{k+1}(-1)^i{n\choose k-i+1}\frac{\ell^i}{i!}.$$ Then $$\frac{(-1)^{n+1}}{n!}\sum_{k=0}^{n-1}A_{n,k}(-1)^k(k+1)!\sum_{i=0}^{k+1}(-1)^i{n\choose k-i+1}\frac{1}{i!}\text{Li}_{-(n+i-k-1)}(1/e).$$ Let $m=i-k$, $$\frac{(-1)^{n+1}}{n!}\sum_{m=-n+1}^1\text{Li}_{-(n+m-1)}(1/e)\sum_{\substack{i-k=m\\0\leq k\leq n-1\\0\leq i\leq k+1}}A_{n,k}(-1)^k(k+1)!(-1)^i{n\choose 1-m}\frac{1}{i!}.$$ This is equivalent to $$\frac{(-1)^{n+1}}{n!}\sum_{m=-n+1}^1\text{Li}_{-(n+m-1)}(1/e)\sum_{i=0}^{n+m-1}A_{n,i-m}(-1)^{i-m}(i-m+1)!(-1)^i{n\choose 1-m}\frac{1}{i!}.$$ Now replace $m$ by $m-n+1$, then $$\frac{1}{n!}\sum_{m=0}^n(-1)^m\text{Li}_{-m}(1/e)\sum_{i=0}^mA_{n,i-m+n-1}(n-m+i)!{n\choose n-m}\frac{1}{i!}.$$ Now the inner sum turns out to be the Stirling numbers of the second kind, $n!S_{m+1}^{(n)}$. Using a closed form for the polylogarithm, $$\sum_{m=0}^n(-1)^mS_{n+1}^{(m+1)}\sum_{k=0}^mk!S_{m+1}^{(k+1)}\frac{1}{(e-1)^{k+1}}.$$ Switching sums, $$\sum_{k=0}^n\frac{k!}{(e-1)^{k+1}}\sum_{m=k}^n(-1)^mS_{n+1}^{(m+1)}S_{m+1}^{(k+1)}.$$ Now using the binomial theorem, $$\frac{(-1)^n}{(e-1)^{n+1}}\sum_{k=0}^nk!\sum_{r=0}^{n-k}{n-k\choose r}e^r(-1)^{k+r}\sum_{m=k}^n(-1)^mS_{n+1}^{(m+1)}S_{m+1}^{(k+1)}.$$ Switching sums again, \begin{eqnarray*} \vartheta_n&=&\sum_{k=0}^\infty\frac{k^n}{k!}B_k\\&=&\sum_{k=0}^n\frac{k!}{(e-1)^{k+1}}\sum_{m=k}^n(-1)^mS_{n+1}^{(m+1)}S_{m+1}^{(k+1)}\\&=&\frac{(-1)^n}{(e-1)^{n+1}}\sum_{r=0}^n\left((-1)^r\sum_{k=0}^{n-r}k!{n-k\choose r}(-1)^k\sum_{m=k}^n(-1)^mS_{n+1}^{(m+1)}S_{m+1}^{(k+1)}\right)e^r, \end{eqnarray*} with $$\frac{e^t}{e^{e^t}-1}=\sum_{k=1}^\infty\frac{t^k}{k!}\vartheta_k.$$
Is there a quicker proof?
Incidentally it's also the case for the Euler numbers that\begin{eqnarray*} \varphi_n&=&\frac{1}{2}\sum_{k=1}^\infty\frac{k^n}{(k-1)!}E_{k-1}(0)\\&=&\sum_{k=0}^n\frac{(-1)^kk!}{(e+1)^{k+1}}\sum_{m=k}^n(-1)^mS_{n+1}^{(m+1)}S_{m+1}^{(k+1)}\\&=&\frac{1}{(e+1)^{n+1}}\sum_{r=0}^n\left(\sum_{k=0}^{n-r}k!{n-k\choose r}(-1)^k\sum_{m=k}^n(-1)^mS_{n+1}^{(m+1)}S_{m+1}^{(k+1)}\right)e^r, \end{eqnarray*} with $$\frac{e^t}{e^{e^t}+1}=\sum_{n=0}^\infty\frac{t^n}{n!}\varphi_n.$$
Let $T_n(x) = \sum_{k=0}^{n} \left\{{ n \atop k }\right\} x^k$ be the Touchard polynomial, where $\left\{{ n \atop k }\right\}$ is the Stirling numbers of the second kind. Then the exponential generating function (EGF) of $T_n(x)$ is given by
$$ \sum_{n=0}^{\infty} \frac{z^n}{n!}T_n(x) = e^{x(e^z-1)}. $$
To better utilized this formula, we define $U_n(x) = \frac{T_{n+1}(x)}{x} = \sum_{k=0}^{n} \left\{{ n + 1 \atop k + 1 }\right\} x^k$. Then the EGF of $U_n(x)$ can be computed as
$$ \sum_{n=0}^{\infty} \frac{z^n}{n!} U_n(x) = \frac{1}{x} \frac{\partial}{\partial z} e^{x(e^z-1)} = e^z e^{x(e^z-1)}. \tag{1} $$
We will also abbreviate $\frac{\partial}{\partial z} = \partial_z$ and likewise for any other partial derivatives. Then
\begin{align*} \frac{e^z}{e^{e^z} - 1} &= \left[ e^z e^{-\partial_w(e^z - 1)} \frac{1}{e^{1 - w} - 1} \right]_{w = 0} \stackrel{(1)}= \sum_{n=0}^{\infty} \frac{z^n}{n!} \left[ U_n(-\partial_w) \frac{1}{e^{1 - w} - 1} \right]_{w = 0}. \end{align*}
In the first step, we utilized the fact that $\partial_w$ is the infinitesimal generator of the translation operator, i.e., $e^{a\partial_w} f(w) = f(w+a)$. (This can also be viewed as the formal way of representing the Taylor expansion.) From this, we obtain:
$$ \bbox[color:navy;padding:5px;border:1px navy dotted;]{\vartheta_n = \left[ U_n(-\partial_w) \frac{1}{e^{1 - w} - 1} \right]_{w = 0} = \sum_{m=0}^{n} (-1)^m \left\{{ n + 1 \atop m + 1 }\right\} \left[ \partial_w^m \frac{e^w}{e - e^w} \right]_{w = 0} } \tag{2} $$
Now we simplify the last sum in $\text{(2)}$. To this end, note that $\frac{e^w}{e - e^w}$ is the EGF of $\left[ \partial_w^m \frac{e^w}{e - e^w} \right]_{w = 0}$. So it suffices to identify the EGF of $\frac{e^w}{e - e^w}$. Indeed,
\begin{align*} \frac{e^w}{e - e^w} = \int_{0}^{\infty} e^w e^{-(e-1)s} e^{(e^w - 1)s} \, \mathrm{d}s \stackrel{(1)}= \sum_{m=0}^{\infty} \frac{w^m}{m!} \int_{0}^{\infty} U_m(s) e^{-(e-1)s} \, \mathrm{d}s \end{align*}
From this, we obtain the integral representation of $\left[ \partial_w^m \frac{e^w}{e - e^w} \right]_{w = 0}$, which then reduces to
\begin{align*} \left[ \partial_w^m \frac{e^w}{e - e^w} \right]_{w = 0} &= \int_{0}^{\infty} U_m(s) e^{-(e-1)s} \, \mathrm{d}s \\ &= \sum_{k=0}^{m} \left\{{ m + 1 \atop k + 1 }\right\} \int_{0}^{\infty} s^k e^{-(e-1)s} \, \mathrm{d}s \\ &= \sum_{k=0}^{m} \left\{{ m + 1 \atop k + 1 }\right\} \frac{k!}{(e - 1)^{k+1}}. \end{align*}
Plugging this back to $\text{(2)}$, we conclude:
$$ \bbox[color:navy;padding:5px;border:1px navy dotted;]{\vartheta_n = \sum_{m=0}^{n} \sum_{k=0}^{m} (-1)^m \left\{{ n + 1 \atop m + 1 }\right\} \left\{{ m + 1 \atop k + 1 }\right\} \frac{k!}{(e - 1)^{k+1}} } \tag{3} $$
This is the same as what OP obtained.