Fermions, different species and (anti-)commutation rules My question is straightforward:
Do fermionic operators associated to different species commute or anticommute? Even if these operators have different quantum numbers? How can one prove this fact in a general QFT?
 A: I think the existing answers to the question are incomplete and confusing, so let me elaborate.
The Hilbert space of two independent fermions of different species (think electron and proton) is the tensor product of the respective spaces, and the field operators for the two fermions act independently on the two factors in the tensor product ($\psi_1\otimes I_2$ and $I_1\otimes \psi_2$ with $I_{1,2}$ being the identity operators on the two small Hilbert spaces). Operators like this commute just by definition, since they concern completely independent degrees of freedom.
This applies to absolutely any quantum fields: as long as they create/annihilate independent degrees of freedom, they commute. Two fermions, two bosons, a fermion and a boson, you name it. The only time anti-commutation comes up is when you write down the relations for one fermion with itself.
However, sometimes you have several "flavors" of fermions that "mix" under the action of some group (think quarks under $SU(3)$). In this case they have to anti-commute for the relations to remain invariant under the action of the mixing group. This is a sign that the better way to think about flavors is as mere components of one larger fermionic field with an extra index, which now doesn't just enumerate the fermions but turns them into a single element of a representation of a group.
To summarize, the "species" of any two fields commute, and the "flavors" of one fermionic field anti-commute (because they are really just components of one bigger fermion).
A: I want to give another perspective. Bosons vs. fermions are usually introduced by considering the many-body wavefunction
$$|\psi(x_1,\ldots,x_N)\rangle$$
in the position basis, where $x_1,\ldots,x_N$ are the positions of $N$ particles. If these $N$ particles are distinguishable, then permuting the $x_j$ should result in a physically equivalent wavefunction, which differs from the original one by a phase. These phases have to be consistently defined, meaning they give a homomorphism from the permutation group $S_N \to U(1)$. It turns out there are two such homomorphisms, a trivial one and one which is $-1$ to the number of swaps. These correspond to bosonic or fermionic particles, respectively.
Now we consider a Hilbert space that includes $N$-particle states for all $N \ge 0$, with a state $|0\rangle$ we call the vacuum. We define the creation operator $a^\dagger(x)$ which creates a particle at $x$ (all particles are indistinguishable still). Then we define
$$|\psi(x_1,\ldots,x_N)\rangle = a^\dagger(x_N)\cdots a^\dagger(x_1)|0\rangle.$$
We see that our permutation action is thus equivalent to commutativity or anticommutativity of the creation operators in the bosonic or fermionic case, respectively.
This permutation symmetry imposes real constraints on the allowed $|\psi(x_1,\ldots,x_N)\rangle$, such as the Pauli exclusion principle for fermions.
However, permutation symmetry is not as impressive with multiple particle types. Suppose we had another species of particles at positions $y_1,\ldots,y_M$ and a joint wavefunction
$$|\psi(x_1,\ldots,x_N,y_1,\ldots,y_M)\rangle.$$
This wavefunction is fixed by $S_N \times S_M$, up to phases, which determine whether each species is bosonic or fermionic. There is no way to swap an $x$ particle with a $y$ particle, however.
We can likewise define another set of creation operators $b^\dagger(y)$ which create a $y$ particle at position $y$. And for which
$$|\psi(x_1,\ldots,x_N,y_1,\ldots,y_M)\rangle = b^\dagger(y_M) \cdots b^\dagger(y_1)a^\dagger(x_N)\cdots a^\dagger(x_1) |0\rangle.$$
We see that the permutation representations of $S_N$ and $S_M$ determine the commutation relations of the two species of creation operators with themselves, but do they determine the commutation relations of them with each other? Actually contrary to several other answers they do not.
To see why, let's study two distinct species of fermion. Distinct physically means that the individual particle numbers, which are schematically
$$N = \int dx a^\dagger(x) a(x) \qquad M = \int dx b^\dagger(x) b(x)$$
are both conserved. This is key to the argument. Note these are Hemitian.
Suppose we now do the usual quantization, where $a$ and $b$ anti-commute. Let us define
$$b' = (-1)^N b.$$
We observe that the $b'$ anti-commute with each other but commute with the $a$'s. Furthermore, we can define the many-body states using the $b'$'s and they only differ by what we had above by a sign. Finally, the time evolution of the $b'$ is equivalent to the time evolution of $b$ because $N$ is conserved.
The reason why some other answers got it wrong thinking about fancy things like $\mathbb{Z}_2$-graded algebras and super-Poisson brackets is that when there are multiple species, there are multiple gradings: in this case a $\mathbb{Z}_2 \times \mathbb{Z}_2$ grading.
In a geometric point of view, which is relevant for bosonization, we would say that neutral fermions all couple to the same spin structure. However, when there are global symmetries, then we can attach charge operators to our fermions for these global symmetries (as we did above), effectively creating multiple spin structures seen by different fermion species.
For more, see this paper: https://arxiv.org/abs/1312.0831
A: Fermionic creation and annihilation operators always satisfy commutation relations with bosonic (or more generally even) operators and anticommutation relations with the fermionic creation and annihilation operators (or more generally odd operators).  This follows from the properties of super Poisson brackets. See Poisson Superalgebra
In particular, the creation operators for distinct, orthogonal modes always anticommute. 
A: Fermions $f_i,\,f_j$ with respective momenta $\pi_i,\,\pi_j$ satisfy the equal-time canonical anticommutation relations $$\left\{\ f_i,\,f_j \right\} = \left\{\ \pi_i,\,\pi_j \right\} = 0,\,\left\{\ f_i\left(t,\,\mathbf{x}\right),\,\pi_j \left(t,\,\mathbf{x'}\right)\right\} = i\hbar \delta_{ij} \delta \left(\mathbf{x},\,\mathbf{x'}\right),$$where the second $\delta$ is a Dirac delta. The $i=j$ special case is a generalisation of a theory of a single fermon $f$ of momentum $\pi$. Why do the $i\neq j$ cases use anticommutators instead of commutators? Because we want our rules to be invariant under $f_i \to \sum_j M_{ij} f_j,\,\pi_i \to \sum_j \left(M^{-1}\right)_{ji} \pi_j$ for invertible choices of the matrix $M$. There's no consistent way to achieve this by using commutators sometimes. A similar explanation is available in terms of the ladder operators.
