Linear functional equations can be solved with Fourier transform.
Let $\lambda_k$ be the roots of the equation $e^\lambda-1=\lambda$. There are infinitely many
such roots. Then

$$f(z)=\sum_k a_ke^{\lambda_k z}$$

is a solution.
Here the sum can be finite, and $a_k$ arbitrary, or the sum can be infinite, and $a_k$ tend to
zero with such speed that the series converges.

The trivial solution is covered, if you interpret the answer correctly.
The equation $e^\lambda-1=\lambda$ has a DOUBLE root at $0$. For the case of a double root
one includes not only the term $e^{\lambda z}$ but also the term $ze^{\lambda z}$.
Similarly for multiple roots $z^ke^{\lambda z}$. So the root $\lambda=0$ exactly covers
the solution $az+b$.

EDIT: Actually all entire solutions can be represented in this way, for the proof I refer to
Gelfond, Calculus of finite differences, MR0342890, Chap. 5 Sect 7, Thm II and Corollaries.
He gives the equation $f'(z)=f(z-1)$ as an example, but your equation is treated similarly.

EDIT2: The problem can be generalized as follows: Let $\omega$ be a distribution with bounded
support. Equation $f\star w=0$, where $\star$ is the convolution, is called a *convolution equation*, and its solutions are called *mean-periodic* functions (fonctions moyenne-periodiques). The case we are discussing is $\omega(x)=\delta(x+1)-\delta(x)-\delta'(x)$
where $\delta$ is the delta function. The theory of mean periodic functions was created by
Delsarte and Schwartz in 1940-s. Schwartz has a general result that all mean-periodic functions
can be obtained as limit of exponential sums, like in this problem.