$\newcommand\R{\mathbb R}\newcommand\C{\mathbb C}$You were almost there.
Indeed, suppose that
$$f(x)=\sum_{j=0}^\infty a_j e^{-jx}\in\R \quad \forall x\in(0,\infty). \tag{1}\label{1}$$
Then
$$g_f(s):=f(-\ln s)=\sum_{j=0}^\infty a_j s^j\in\R \quad \forall s\in(0,1).$$
So, the radius of convergence of the power series $\sum_{j=0}^\infty a_j z^j$ is $1$. So, $g_f$ can be extended to the analytic function $G_f$ on the open unit disc
$D$ in $\C$ given by the formula $G_f(z)=\sum_{j=0}^\infty a_j z^j$ for $z\in D$.
Vice versa, if the function $g_f$ can be extended to an analytic function $G_f$ on
$D$, then $G_f(z)=\sum_{j=0}^\infty a_j z^j$ with $a_j=G_f^{(j)}(0)/j!$ and hence \eqref{1} holds.
Thus, \eqref{1} holds for some $a_j$'s if and only if the function $g_f$ can be extended to an analytic function on $D$.
Similarly, \eqref{1} holds for some $a_j$'s if and only if the function $g_f$ can be extended to a real-analytic function on the interval $(-1,1)$.
For an illustration, suppose that $f(x)=\dfrac1{1-e^{-x}}$ for $x\in(0,\infty)$. Then $g_f(s)=f(-\ln s)=\dfrac1{1-s}$ for $s\in(0,1)$. Clearly, the function $g_f$ can be extended to the analytic function $G_f$ on $D$ given by the formula
$$G_f(z):=\dfrac1{1-z}=\sum_{j=0}^\infty a_j z^j=\sum_{j=0}^\infty z^j$$
with $a_j=G_f^{(j)}(0)/j!=1$. So,
$$f(x)=\sum_{j=0}^\infty a_j e^{-jx}=\sum_{j=0}^\infty e^{-jx}$$
for $x\in(0,\infty)$.
