What is the fundamental reason of the fermion doubling? Recall that the fermion doubling is the problem in taking the $a \to 0$ limit of a naively discretized fermionic theory (defined on a lattice with lattice spacing $a$). After such a limit one finds themselves with an additional amount (precisely $2^d$) of fermionic fields. One can fix this by considering different discretizations of the action that make unwanted fields decouple in the continuum limit. The downside is that the additional terms have to spoil some nice features of the theory (chiral symmetry, locality, lattice symmetry, ...).
Now, I wonder what is the true reason for the appearance of new fields. Is it the fermionic nature of the theory? (In other words, is a similar problem ruled out for bosonic fields?) And do all (naive?) fermionic theories (that is, independent of the continuum form of the action) suffer from this problem?
More generally, how can one tell a priori what will the field content of a lattice theory in the continuum limit be? Or is the field content fundamentally a continuum limit property that has to be calculated?
 A: What is the fundamental reason of the fermion doubling?
Answer: there is no fundamental reason of the fermion doubling.
Adding proper lattice interaction can always get rid of the doublers, and it works for both abelian and non-abelian gauge groups, as long as the resulting chiral theory is free of all anomalies. (See Xiao-Gang Wen arXiv:1305.1045 and Yi-Zhuang You Cenke Xu arXiv:1412.4784).
A: A density matrix view of the problem:
When $p_\mu$ is near zero, one is considering a momentum near zero and the discrete lattice works fine. It's when $p_\mu$ is near $\pm \pi/a$, or more generally, near $n\pi/a$ for $n\neq 0$ an integer, that one finds problems. Instead of having very large momenta, these values of $p_\mu$ essentially give momenta as small as those near $p_\mu=0$, but with signs in the $\gamma$ matrices negated.
This is a type of aliasing problem. As with the usual aliasing, the problem at $\pm \pi/a$ goes away when one makes the sample rate faster, that is, replaces $a$ with a smaller value. And just as with aliasing, making $a$ smaller does not eliminate aliasing entirely, but instead pushes the problem to a higher frequency.
This drawing shows the usual aliasing effect. Note that the high frequency black signal (which corresponds to a high momentum) appears as a low frequency red signal (with low momentum):
 

The difference with the usual aliasing is what happens when $p_\mu = (2n+1)\pi/a$. These are the values that give a continuous Dirac equation with negated gamma matrices. To understand better what is going on here, let's consider the density matrix form. Density matrices avoid unphysical complex phases. I'll work in 3+1 dimensions.
One obtains a density matrix $\rho$ by multiplying a ket by a bra:
$$\rho = |a\rangle\langle a|.$$
The reason there are four spin-1/2 particles in 3+1 dimensions is that there are four primitive (i.e. trace=1) solutions to the idempotency equation:
$$\rho^2 = \rho$$
One has latitude in how one chooses the four states. In general, one chooses two elements of the Dirac algebra that (a) square to unity, (b) commute, and (c) are independent. These are called a "complete set of commuting roots of unity."
The commuting roots of unity are operators; one chooses them according to what operators one wishes to diagonalize. Following the wikipedia article on the construction of Dirac spinors, if we choose z-spin and charge $Q$, our complete set of commuting roots of unity is:
$$\sigma_z = i\gamma^1\gamma^2,\;\;\; Q = -\gamma^0$$
The four independent states are then:
$$\rho = (1\pm \sigma_z)(1\pm Q)/4.$$
To get the spinors from a density matrix, one chooses a nonzero column and normalizes. Thus spinors and density matrices are alternative mathematical representations of wave functions; neither is more fundamental.

If we discretize the density matrix, we will end up with the usual aliasing problem. From the point of view of lattice type calculations this is acceptable; there will be no duplicated particles. But spinors carry an extra degree of freedom; the arbitrary complex phase. This makes their aliasing behavior more complicated.
So consider what happens for a frequency a little larger than $\pi/2$. In the following illustration we color every other sample red or blue:

In the above, the density matrix will see this frequency appropriately as a high frequency. But with a spinor, we have arbitrary complex phase freedom. So we can negate the blue dots; the result is a low frequency. Thus the arbitrary complex phases of spinors naturally give aliasing problems at half frequencies.
A: The fermion doubling is manifested through the existence of extra poles in the Dirac propagator on the lattice. These poles cannot be made to disappear at the continuum limit. (The number of doublers can be reduced by different discretizations but not eliminated at all, this is essentially the Nielsen-Ninomiya theorem).
The reason for the fermion doubling lies in the existence of chiral anomaly. This anomaly exists in the continuum limit due to the chiral non-invariance of the path integral measure and not because of the non-invariance of Lagrangian. In the lattice formulation of a chiral theory based on a discretization of the Lagrangian, the anomaly is absent and the lattice formulation generates the extra species just to cancel this anomaly in the continuum limit. Since the axial anomaly exists in nature $\pi^0 \rightarrow \gamma \gamma$, this situation is unacceptable.
The fermion doubling problem is an artifact of the realization of the theory by means of quarks where the axial anomaly is not present in the lagrangian but in the path integral measure. There are approaches of other types of discretizations such as by means of fuzzy spaces in noncommutative geometry, where the quarks are not the basic fields of the theory. In these approaches, the fermion doubling problems do not exist.
Update
This is an update referring to Marek's first comment.
For fermions, the axial anomaly can be recovored only by a non-trivial regularization of the path integral measure of the fermions. 
Any finite dimensional approximation of this measure as a product of Berezin-Grassmann integrals does not produce the axial anomaly.
This is the reason why the Lattice regularization does not produce the axial anomaly and as a consequence the phenomenon of doubling occurs where the
different species of doublers are of opposite chirality to cancel the anomaly. This is a property of fermion fields represented by Grassmann variables.
In effective field theories (where the basic fields are pions) such as sigma models, the axial anomaly is manifested through a Wess-Zumino-Witten term in the Lagrangian, but these theories are not even 
perturbatively renormalizable and I think this is the reason why they were not put on a lattice.
One approach that I know of which solves the fermion doubling problem is the regularization by means of a fuzzy space
approximation of the space-time manifold. The philosophy of this approach is explained in the introduction of Mark's Rieffel's article. The resolution of the fermion doubling using this approach is given nicely in Badis Ydri's thesis. (There is also more recent work on the subject by A. Balachandran and B. Ydri in the arxiv).
The main idea is that the Poisson algebra of functions over certain spaces (such as the two sphere, and more generally coadjoint orbits of compacl Lie groups) can be approximated by finite dimensional matrices. These approximations are called fuzzy spaces. Gauge fields and Fernmions can be constructed on these fuzzy spaces, which have the correct continuum limit when the matrix dimensions become infinite. This formulation contains the axial anomaly inherently, thus it is free of the doubling problem.
The only drawback that I can see of this approach is that it is applicable only to some special 4-dimensional manifolds such as $CP^2$ or $S^2 \times S^2$, because the fuzzy manifold is required to be Berezin-quantizable.
