This is currently a partial answer, refining the idea given by @reuns.
The series $\Phi(s)=\sum_{n=1}^\infty\ e^{-n^s}$ converges iff $s>0$ is real. Using the Cahen–Mellin integral $$e^{-x}=\frac{1}{2\pi i}\int_{c-i\infty}^{c+i\infty}\Gamma(z)x^{-z}\,dz\qquad(x,c>0)$$ with $x=n^s$ and $c>1/s$, we get $$\Phi(s)=\frac{1}{2\pi i}\int_{c-i\infty}^{c+i\infty}\Gamma(z)\zeta(sz)\,dz.$$
For $0<s<1$, the integrand tends to $0$ rapidly enough when $z\to\infty$ in the half-plane $\Re z\leqslant c$ and out of a neighborhood of the line $L=\{z : \Im z=0\wedge\Re z\leqslant 1/s\}$. This allows us to deform the path of integration, making it encircle $L$, and we see that $\Phi(s)$ is equal to the (infinite) sum of residues of the integrand at its poles (which are $z=1/s$ and $z=-n$ for nonnegative integers $n$). Computing these, we get $$\Phi(s)=\Gamma\left(1+\frac1s\right)+\sum_{n=0}^\infty\frac{(-1)^n}{n!}\zeta(-ns).$$
This series converges for complex $s\neq 0$ with $\Re s<1$ (at least; the singularities at $s=-1/n$ for $n\in\mathbb{Z}_{>0}$ are removable), and gives the analytic continuation of $\Phi(s)$ in this region.
The remaining question is whether we can extend this region further.