This is a partial answer. I will show that for $a>0$ small enough
$$
f(x) = x(1+x^n)^n - a(1+x^n)^n + a^nx
$$
will have exactly one positive real zero and for $a>0$ large enough $f$ will have exactly three positive real zeros.
The case when $a$ is small.
When $a = 0$ the zeros of $f$ are the $n^\text{th}$ roots of $-1$, each with multiplicity $n$, and $x=0$, with multiplicity $1$.
Figure 1: The zeros of $f$ for $n=5$ when $a=0$. The unit circle is shown in blue for reference.
Zeros of polynomials are continuous functions of their coefficients, which means that for any $\epsilon > 0$ and any $n^\text{th}$ root of $-1$, say $\zeta$, $f$ will have exactly $n$ zeros in the disk $|x-\zeta| < \epsilon$ for $a$ sufficiently small. In particular, none of the $n^2$ zeros of $f$ clustered around the $n^\text{th}$ roots of $-1$ will lie on the real axis for $a$ sufficiently small.
Figure 2: The zeros of $f$ for $n=5$ with $0 < a < 0.9$. As $a$ increases from $0$, exactly $5$ zeros erupt from each of the $5^\text{th}$ roots of $-1$. We can also see the zero which was previously located at $x=0$ move to the right. Also, two of the zeros coming from the primitive root of $-1$ and its conjugate travel toward the real axis.
Now $f(0) = -a$ and $f(x) \to +\infty$ as $x \to +\infty$, so $f$ will have at least one positive real zero when $a > 0$. By the above remark this positive real zero will be unique for $a > 0$ sufficiently small.
Aside 1: The implicit function theorem can be invoked to see that this zero, which we will call $x_0 = x_0(a)$, is actually a differentiable function of the parameter $a$ as long as $a$ is sufficiently small. With this in mind, by differentiating the equation
$$
x_0(1+x_0^n)^n - a(1+x_0^n)^n + a^nx_0 = 0
$$
with respect to $a$ it can be shown that $x_0'(0) = 1$. This implies that $x_0(a) \sim a$ as $a \to 0$.
The case when $a$ is large.
We noted above that $f(0) < 0$ for all $a > 0$. Further, $f(x) \to \infty$ as $a \to \infty$ for fixed $x>0$, which means that $f$ has a positive real zero which tends toward the point $x=0$ as $a \to \infty$.
Aside 2: For $x>0$ and $a>0$ we can show that for any $0 < \epsilon < 1$ we will have $f\left(\epsilon a^{1-n}\right) < 0$ and $f\left(a^{1-n}\right) > 0$ for $a$ sufficiently large provided $n > 1$. It follows that $f$ has a zero in the interval $(\epsilon a^{1-n},a^{1-n})$ for $a>0$ sufficiently large. If we call this zero $\hat{x}_0(a)$ then this implies that $\hat{x}_0(a) \sim a^{1-n}$ as $a \to \infty$ with $a > 0$.
Numerically we can see that, for $a>0$ large, $f$ has a zero which is approximately equal to $a$ (see Figures 3 and 4). Indeed we see that, for $n>1$,
$$
\begin{align}
\frac{f(ax)}{a^{n^2+1}} &= x\left(a^{-n}+x^n\right)^n - \left(a^{-n} + x^n\right)^n + a^{n-n^2}x \\
&\to x^{n^2+1} - x^{n^2} \\
&= x^{n^2}(x-1) \tag{$*$}
\end{align}
$$
as $a \to \infty$, where the convergence is uniform with respect to $x$ in compact subsets of $\mathbb{C}$. From this we may conclude that $f(ax)$ has a zero arbitrarily close to $1$ for $a$ large enough, and so $f(x)$ has a zero, call it $x_\infty(a)$, for which $x_\infty(a) \sim a$ as $a \to \infty$.
To prove that $x_\infty$ is real we will employ an approach similar to the one we used in Aside 2. Indeed for $a>0$ we have $f(a) = a^{n+1} > 0$ and, for $0 < \epsilon < 1$,
$$
f(\epsilon a) = a (\epsilon - 1) \Bigl(1 + (\epsilon a)^n\Bigr)^n + \epsilon a^{n+1} < 0 \tag{$**$}
$$
for $a > 0$ sufficiently large provided $n > 1$. Thus $f$ has a zero in any interval $(\epsilon a,a)$ for $a > 0$ sufficiently large.
Figure 3: The zeros of $f$ for $n=5$ and $1.2 < a < 3$. In addition to the zeros continuing to radiate outwards from the $5^\text{th}$ roots of $-1$ we also see a zero tending toward the point $x=0$ and a large zero quickly moving to the right on the real axis. These are the zeros asymptotic to $a^{1-n}$ and $a$, respectively.
So far we have found two of the $n^2+1$ zeros of $f$. It turns out that the remaining $n^2-1$ zeros of $f$ are asymptotically evenly spaced on the circle of radius $|x| = a^{1/(n+1)}$ as $a \to \infty$.
We have
$$
\begin{align}
\frac{f\left(a^{1/(n+1)} x\right)}{a^{n/(n+1)-n-1}} &= a^{-n/(n+1)} x \left(a^{-n/(n+1)} + x^n\right)^n - \left(a^{-n/(n+1)} + x^n\right)^n + x \\
&\to x - x^{n^2} \\
&= x\left(1-x^{n^2-1}\right)
\end{align}
$$
as $a \to \infty$, where the convergence is uniform with respect to $x$ in compact subsets of $\mathbb{C}$. It follows that $f\left(a^{1/(n+1)} x\right)$ has a zero arbitrarily close to $\omega$ for $a$ large enough, where $\omega$ is any $(n^2-1)^\text{th}$ root of unity. Consequently $f(x)$ has a zero which is asymptotic to $\omega a^{1/(n+1)}$ as $a \to \infty$ for each $\omega$.
In particular, $f$ has a zero which is asymptotic to $a^{1/(n+1)}$ as $a \to \infty$. To show that this is real we proceed as before, observing that $f\left(a^{1/(n+1)}\right) < 0$ and, for any $0 < \epsilon < 1$, $f\left(\epsilon a^{1/(n+1)}\right) > 0$ for $a > 0$ sufficiently large. Thus $f$ has a zero in the interval $\left(\epsilon a^{1/(n+1)}, a^{1/(n+1)}\right)$ for $a>0$ sufficiently large.
Figure 4: The zeros of the rescaled polynomial $f\left(a^{1/(n+1)} x\right)$ for $n=5$ and $0 < a < 5$. We can see the zeros which originally erupted from the $n^\text{th}$ roots of $-1$ (prior to rescaling) travel toward points on the unit circle. This rescaled polynomial also has a large real zero which will tend to $\infty$ as $a \to \infty$.
Figure 5: The last shown state of the zeros displayed in Figure 4. Specifically, these are the zeros of the rescaled polynomial $f\left(a^{1/(n+1)} x\right)$ for $n=5$ and $a = 5$.
Aside 3: I guessed the scale factor $a^{1/(n+1)}$ based on numerical experiments but in principle it could have been deduced by trying to balance the two largest terms of $f(a^q x)$. We showed in $(*)$ that we would need to take $q < 1$, and by noticing that
$$
\frac{f(x)}{a^n} \to x
$$
as $a \to \infty$ provided $n>1$ we would also know that we must take $q > 0$.
From the discussion above we may deduce that the remaining $n^2-1$ zeros do not lie on the real axis for $a>0$ sufficiently large. Thus $f$ has exactly three positive real zeros for $a>0$ sufficiently large.
Summary of the results.
We will suppose always that $a>0$.
For $a$ small the polynomial $f$ has exactly one real zero which is asymptotic to $a$ as $a \to 0$. The remaining $n^2$ zeros tend to the $n^\text{th}$ roots of $-1$ as $a \to 0$.
For $a$ large the polynomial has exactly three real zeros, one asymptotic to $a^{1-n}$, another asymptotic to $a$, and the last asymptotic to $a^{1/(n+1)}$ as $a \to \infty$. The remaining $n^2-1$ zeros are asymptotic to $\omega a^{1/(n+1)},$ where $\omega$ is an $(n^2-1)^\text{th}$ root of unity.
Let $f(x)$ be any quartic polynomial with coefficients from $\{ -1, +1 \}$. Replacing $f(x)$ by $-f(x)$ if necessary, we can assume $f(x)$ is monic. i.e.
$$f(x) = x^4 + ax^3 + bx^2 + cx + d\quad\text{ with }\quad a,b,c,d \in \{ -1, +1 \}$$
If $f(x)$ has $4$ real roots $\lambda_1,\lambda_2,\lambda_3,\lambda_4$, then by Vieta's formula, we have
$$\sum_{i=1}^4 \lambda_i = -a, \sum_{1\le i < j\le 4} \lambda_i\lambda_j = b
\quad\text{ and }\quad\prod_{i=1}^4 \lambda_i = d$$
Notice
$$\sum_{i=1}^4 \lambda_i^2 = \left(\sum_{i=1}^4\lambda_i\right)^2 - 2\sum_{1\le i < j \le 4}\lambda_i\lambda_j = a^2 - 2b = 1 -2b$$
Since $\sum_{i=1}^4 \lambda_i^2 \ge 0$, we need $b = -1$. As a result, $$\sum_{i=1}^4 \lambda_i^2 = 3$$
By AM $\ge$ GM, this leads to
$$\frac34 = \frac14\sum_{i=1}^4 \lambda_i^2 \ge \left(\prod_{i=1}^4 \lambda_i^2\right)^{1/4} = (d^2)^{1/4} = 1$$
This is impossible and hence $f(x)$ cannot has 4 real roots.
Best Answer
Let $E_n(x):=\sum_{k=0}^n\,\frac{x^k}{k!}$ for $n=0,1,2,\ldots$. We shall prove that $E_n(x)$ has no real roots if $n$ is even, and $E_n(x)$ has exactly one real root, which is simple, if $n$ is odd.
Suppose that $n$ is even. Clearly, $E_n(x)$ has no roots in $\mathbb{R}_{\geq 0}$. By Taylor's Theorem, we have $\exp(x)=E_n(x)+R_n(x)$, where the remainder term is given by $$R_n(x)=\int_0^x\,\frac{\exp^{(n+1)}(t)}{n!}\,(x-t)^n\,\text{d}t=\int_0^x\,\frac{\exp(t)}{n!}\,(x-t)^n\,\text{d}t\,.$$ If $x<0$, then $$R_n(x)=-\int_0^{|x|}\,\frac{\exp(-t)}{n!}\,|x+t|^n\,\text{d}t<0\,.$$ That is, $$E_n(x)=\exp(x)-R_n(x)>\exp(x)>0$$ for all $x<0$. That is, $E_n(x)$ has no negative roots either; i.e., $E_n(x)$ has no real roots.
If $n$ is odd, then $E'_n(x)=E_{n-1}(x)$ has no real roots. Thus, $E_n(x)$ can have at most one real root, due to Rolle's Theorem. Clearly, $E_n(x)$ has a real root, being a polynomial in $\mathbb{R}[x]$ of an odd degree. Consequently, $E_n(x)$ has exactly one real root, which is simple.