The one direction is clear, if $f = g^n$, and $z_0$ is a zero of $f$, then $z_0$ is a zero of $g$, say of order $k$, so in a neighbourhood of $z_0$ we have $g(z) = (z-z_0)^k\cdot h(z)$ with a holomorphic nonzero $h$. Then in that neighbourhood $f(z) = (z-z_0)^{nk}\cdot h(z)^n$, and of course $h(z)^n$ is nonzero there. So indeed $n$ divides the order of the zero of $f$ at $z_0$.
Conversely, let all zeros of $f$ have orders divisible by $n$. The case $f \equiv 0$ is trivial, so in the following we assume $f\not\equiv 0$. Consider the entire meromorphic function
$$h_0(z) = \frac{f'(z)}{f(z)}.$$
Since $f$ is entire, the poles of $g$ are exactly the zeros of $f$. Let $P$ be the pole set of $g$ and $U = \mathbb{C}\setminus P$.
All poles of $h_0$ are simple, with the residue of $h_0$ in $z_0$ being the order of the zero of $f$ at $z_0$: If $f(z) = (z-z_0)^m\cdot k(z)$ with a holomorphic nonzero $k$ in a neighbourhood of $z_0$, then
$$\frac{f'(z)}{f(z)} = \frac{m(z-z_0)^{m-1}\cdot k(z) + (z-z_0)^m\cdot k'(z)}{(z-z_0)^m\cdot k(z)} = \frac{m}{z-z_0} + \frac{k'(z)}{k(z)}$$
in that neighbourhood, and $\frac{k'}{k}$ is holomorphic there since $k$ is nonzero.
All residues of $h_0$ are multiples of $n$, hence the entire meromorphic function $h = \frac{1}{n}\cdot h_0$ has integer residues.
Now choose an $a \in U$, and for every $z\in U$ choose a - piecewise continuously differentiable - path $\gamma_z \colon [0,1] \to U$ with $\gamma_z(0) = a$ and $\gamma_z(1) = z$. Then define $\lambda \colon U \to \mathbb{C}$ by
$$\lambda(z) := \int_{\gamma_z} h(\zeta)\,d\zeta.$$
The function $\lambda$ depends on the choice of the path $\gamma_z$, but if we choose a different path $\delta_z$ from $a$ to $z$, then the composition $\gamma_z\delta_z^{-1}$ is a closed path in $U$, and $n(\gamma_z\delta_z^{-1},\omega) \neq 0$ for only finitely many $\omega\in P$. Hence
\begin{align}
\int_{\gamma_z} h(\zeta)\,d\zeta - \int_{\delta_z} h(\zeta)\,d\zeta
&= \int_{\gamma_z\delta_z^{-1}} h(\zeta)\,d\zeta\\
&= 2\pi i \sum_{\omega\in P} n(\gamma_z\delta_z^{-1},\omega)\operatorname{Res}(h;\omega)\\
&\in 2\pi i\cdot \mathbb{Z},
\end{align}
so a different choice of paths alters $\lambda(z)$ by an integer multiple of $2\pi i$, whence $g_0 \colon U\to \mathbb{C}$ given by
$$g_0(z) = \exp(\lambda(z))$$
is independent of the choice of paths. The independence of the choice of paths shows that $g_0$ is holomorphic: For $w\in U$ and $r > 0$ such that $D_r(w) \subset U$, we can define
$$\tilde{\lambda}(z) = \lambda(w) + \int_w^z h(\zeta)\,d\zeta$$
on $D_r(w)$, which is holomorphic, and we have
$$\exp(\tilde{\lambda}(z)) \equiv \exp(\lambda(z))$$
on $D_r(w)$ by the independence of $g_0$ from the choice of paths. Also, this shows that $g_0'(z) = \frac{f'(z)}{nf(z)}\cdot g_0(z)$.
Now we see that $f\cdot g_0^{-n}$ is constant on $U$ by differentiation:
$$(f\cdot g_0^{-n})'(z) = f'(z)\cdot g_0^n(z) - nf(z)\cdot g_0^{n-1}(z)\cdot g_0'(z) = f(z)g_0^n(z)\biggl( \frac{f'(z)}{f(z)} - n\frac{g_0'(z)}{g_0(z)}\biggr) = 0.$$
Hence, for a suitable $c\in \mathbb{C}\setminus \{0\}$, $g(z) := c\cdot g_0(z)$ is a holomorphic $n$-th root of $f$ on $U$.
But for $\omega\in P$ we have
$$\lim_{z\to \omega} g(z) = 0,$$
hence these are all removable singularities and $g$ extends to an entire function.
Is an argument of the same form possible using a Weierstrass product in place of the polynomial here? I'd hate to bring in such a high-tech object for this sort of problem, but entire functions can be weird so perhaps this is necessary.
– Dorebell Aug 12 '15 at 23:43