Updated I think that for non-polynomials we need to restrict to non-negative $x$ (or at least $x \gt -1$). With this restriction, are there discrete analytic functions which are not natural functions? I think not, at least under rather lax conditions.
Consider an arbitrary expansion $f(x)=\sum_0^{\infty}a_k \binom{x}{k}$ with the $a_k$ real. It is defined for all non-negative integer $x=n$ since only the first $n+1$ terms are non-zero. However for negative integral $x$ we have $\binom{x}{k}=(-1)^k\binom{|x|+k-1}{k}$ so divergence is quite possible for negative $x$ values. I think that as long as $\lim_{k \to \infty}\frac{a_k}{k!}=0$ (or at least if $a_k$ has at worst polynomial growth) then we also have $f(x)$ convergent for all real $x \ge 0$ and, by your definition, $f(x)$ is discrete-analytic (on that range). Then $f$ is determined by it's values at the non-negative integers.
Here is a slightly more formal version of my previous comments (which only deal with the first sentence of your post, which isn't even the main question.) I like it well although it is not original. It has that nice quality of seeming mystifying until the moment when it seems trivially obvious.
Consider a member $\mathbf{y}=(y_0,y_1,y_2,\cdots)$ of the space of sequences. We also write $\mathbf{y}(x)=y_x$ for non-negative integer $x$. Aside from the identity operator $I$ we have the difference operator with $\Delta\mathbf{y}=(y_1-y_0,y_2-y_1,\cdots)$ and the shift operator with $S\mathbf{y}=(y_1,y_2,y_3,\cdots)$. So $y_n=\left(S^n\mathbf{y}\right)(0).$
Since $S=I+\Delta,$ and the operators commute, we have $S^n=(I+\Delta)^n=\sum_{k=0}^n\binom{n}{k}\Delta^k$ and, as you said, $$y_n=\sum_{k=0}^{\infty}\binom{n}{k}\left(\Delta^k\mathbf{y}\right)(0)=\sum_{k=0}^{\infty}a_k\binom{n}{k}$$ for the coefficients $a_k=\left(\Delta^k\mathbf{y}\right)(0).$ Also, $$\mathbf{y}=\sum_{k=0}^{\infty}a_k\mathbf{\binom{x}{k}}$$ Where, for example, $\mathbf{\binom{x}{3}}=(0,0,0,1,4,10,20,\cdots)$ is a sequence.
Note too that the $\mathbf{\binom{x}{k}}$ are a natural basis for the space of sequences with $\Delta\mathbf{\binom{x}{k}}=\mathbf{\binom{x}{k-1}}$ and $S\mathbf{\binom{x}{k}}=\mathbf{\binom{x}{k}}+\mathbf{\binom{x}{k-1}}$ for $k \ge 1,$ while $\mathbf{\binom{x}{0}}$ is sent to the $0$ sequence by $\Delta$ and itself by $S$. Accordingly the shift of $\mathbf{y}$ is $$S\mathbf{y}=\sum_{k=0}^{\infty}(a_k+a_{k+1})\mathbf{\binom{x}{k}}. $$
Going back to the expansion for $\mathbf{y}$, we then have a corresponding discrete-analytic function $$f(x)=\sum_{k=0}^{\infty}a_k\binom{x}{k}$$ defined for all non-negative real $x.$ Here convergence is not an issue (under mild conditions) and $\mathbf{y}$ is the restriction of $f$ to the integers. Also, the discrete-analytic function corresponding to $S\mathbf{y}$ is $$\sum_{k=0}^{\infty}(a_k+a_{k+1})\binom{x}{k}. \tag{*}$$
For $f(x)$ to satisfy the definition of a natural function we need the shift $f(x+n)$ to be discrete-analytic for all positive $n.$ It is sufficient to establish this just for $n=1$ provided that we do this for all discrete-analytic functions.
But $$f(x+1)=\sum_{k=0}^{\infty}a_k\binom{x+1}{k}=a_0\binom{x+1}{0}+\sum_{k=1}^{\infty}a_k\left( \binom{x}{k-1}+\binom{x}{k}\right).$$ Thus $$f(x+1)=\sum_{k=0}^{\infty}(a_k+a_{k+1})\binom{x}{k}. \tag{**}$$
To show that $f(x+1)$ is indeed discrete-analytic, examine the sequence arising from the restriction of $g(x)=f(x+1)$ to the non-negative integers. This restriction is $S\mathbf{y}$ so, comparing $( * )$ and $( ** )$, we are done.