The Majorana fermions can be written in terms of usual (spinless) fermions by the relation $ f_j = c_j^A + i c_j^B $, where $A$ and $B$ refer to sublattice and $j$ to the unit cell. This leads to a spinless superconducting Bogoliubov-deGennes (BdG) Hamiltonian, which means it includes not only fermionic hopping terms, but also terms that create and destroy local fermion pairs. This is most easily written in $k$-space.
\begin{align}
H &= \sum_{\left<ij\right>^a} J^a i c_i c_j
= \sum_{\mu} i \left[ J^z c_\mu^A c^B_\mu + J_x c_\mu^A c^B_{\mu+n_1} + J_y c_\mu^A c^B_{\mu+n_2} \right] \\
&= \sum_k \left[f^\dagger_k f_k 2 \text{Re} \Gamma_k - J_z + \left(f_k f_{-k} - f^\dagger_{-k} f^\dagger_k \right) i \text{Im} \Gamma_k \right],
\end{align}
where $\Gamma_k = J^z + J^x e^{i \vec{n}_1 \cdot \vec{k}} + J^y e^{i \vec{n}_2 \cdot \vec{k}}$.
Note that the Hamiltonian does not take the form that you wrote above with $i$ and $j$ lattice sites. The superconducting terms appear on the x and y bonds due to the fact that we chose to pair the Majoranas along the z-bonds (the two sites in a single unit cell were chosen to be connected by a z-bond). These fermions' number is not conserved (although they are modulo 2). However, if we now diagonalize this Hamiltonian by a Bogoliubov transformation we get the following form.
$$H = \sum_k \epsilon_k \left( a^\dagger_k a_k - \frac{1}{2} \right),$$
where $\epsilon_k = 2|\Gamma_k|$. These are quasiparticles, and the excitations of the system are defined by creating and destroying these. The Hamiltonian does conserve the quasiparticle number $\mathcal{N} = \sum_k a_k^\dagger a_k$. In the ground state the number of quasiparticles is zero, and in the highest energy state in the zero-flux Fock space the number is equal to the number of lattice sites since there is an excitation for every $k$ value.
The relation to the original spin model is a difficult one due to the non-local nature of the quasiparticles.
Response to Hamed's comment: I think some comments like "double spectrum" and "taking only the positive eigenvalue" come from the analogy with diagonalizing matrices. The first step to "diagonalize" the Hamiltonian into the form I wrote above is is to write the Hamiltonian in Nambu matrix form. The standard basis is $(f_k,f^\dagger_{-k})^T$, but the the Majorana basis $(c_k^A,c_k^B)$ chosen used Kitaev can save you a step in this case. Then the remaining step is just standard unitary matrix diagonalization due to the commutation relations of the resulting fermions. The resulting matrix is $$H = \frac{1}{2} \sum_k\left( \begin{array}{c} a_k \\ a^\dagger_{-k} \end{array} \right)^\dagger\left( \begin{array}{cc} {\epsilon_k} & 0 \\ 0 & - {\epsilon_k} \end{array} \right) \left( \begin{array}{c} a_k \\ a^\dagger_{-k} \end{array} \right), $$
which reproduces the previous expression for the Hamiltonian due to the commutation relations. In this sense there are two eigenvalues, and the true energy of the quasiparticles corresponds to their difference multiplied by the overall factor of one-half, or the higher energy eigenvalue (due to the cleverly chosen one-half, following Kitaev). In short the spectrum is repeated because the BdG form splits the Hamiltonian into separate $k$ and $-k$ parts and the symmetry of $\epsilon_k$ under $k\to -k$ makes these two parts redundant. This is special to the case of a BdG Hamiltonian with particle-hole symmetry, and contrasts with the case of electrons filling a band where each eigenvalue of the matrix represents a different excitation.
