Why fermions have a first order (Dirac) equation and bosons a second order one? Is there a deep reason for a fermion to have a first order equation in the derivative while the bosons have a second order one? Does this imply deep theoretical differences (like space phase dimesion etc)?
I understand that for a fermion, with half integer spin, you can form another Lorentz invariant using the gamma matrices $\gamma^\nu\partial_\nu $, which contracted with a partial derivative are kind of the square root of the D'Alembertian $\partial^\nu\partial_\nu$.  Why can't we do the same for a boson? 
Finally, how is this treated in a Supersymmetric theory? Do a particle and its superpartner share a same order equation or not?
 A: You have hit on the exact superficiality that drives people into supersymmetry theories. There is nothing fundamentally absolute about a 2-ODE governing the dynamical evolution of a ``particle''. Square your Dirac equation, and you get a 2-ODE (KG), but the former is still more illuminating and contains as much information as the latter. In fact, the requirement of spin and dimensionality (N=4, simplest realizable case) follows naturally when you adopt a critical perspective towards a 2-ODE (it gave Dirac, probability densities that weren't positive definite and that made no physical sense, unless you start inventing new conjectures (e.g. hole theory)). (You know these stories, don't you. Otherwise read Bjorken, Drell - that won't be too many details, but just enough).
let us take the simplest N=4 case (that's the lowest dimension in which you can realize that group structure. But of course, higher ones are possible.) We can make `physical' sense of this dimensionality by recognizing that the four components of the bi-spinor carry information regarding (spin up, spin down)*(positive energy, negative energy) parts. But again that dimensionality was forced by the group structure and is NOT a matter of any fundamental distinction between a particle with spin and a spinless particle. Square the operator part of the Dirac equation, and you will get a 2-ODE for each of these components.
So, this entire formulation creates this misleading impression that dimensionality of space creates any fundamental distinction between fermions and bosons. How this gets treated in SUSY is a tad too detailed to be mentioned here (plus, it will kill the fun, when you eventually learn about it from reading a textbook). I genuinely recommend W. Seigel's SUSY book and would only leave you with a teaser - not just that there is no absolute distinction b/w fermions and bosons, the spectrum is much richer; ever heard of `anyons'? That's OK, you'll find out in due course of time. Happy leaning :)
PS - Like every conscientious physicist, I leave you with a word of WARNING. SUSY is a beautiful theory, could be a great step in the right direction if it is proved correct. BUT DON'T TREAT IT AS A BIBLE, UNTIL THERE IS ANY EXPERIMENTAL EVIDENCE IN FAVOR OF IT, no matter how beautiful it may be.     
A: Let me note that the Dirac equation is a system of first order partial differential equations (PDE) for 4 components of the spinor. You can present the Klein-Gordon equation as a system of first order PDEs for 5 unknowns ($\psi$ and its partial derivatives with respect to $t,x,y,z$: $u=\psi_{,t}$,$v=\psi_{,x}$,$w=\psi_{,y}$, $r=\psi_{,z}$. On the other hand, the Dirac equation is generally equivalent to one fourth order equation for one unknown (see my article http://akhmeteli.org/wp-content/uploads/2011/08/JMAPAQ528082303_1.pdf (JOURNAL OF MATHEMATICAL PHYSICS 52, 082303 (2011))
A: Let me give an answer i find very good and based on a small textbook on relativistic quantum mechanics (Relativistic Quantum Mechanics, S Trahanas, in greek)
First one can start with the known data of particles with integral spin (refered as bosons, obeying Bose-Einstein statistics) and particles with half-integral spin (refered as fermions obeying Fermi-Dirac statistics, or Pauli's exclusion principle)
One then can try to find a relativistic quantum equation (effectively a Schrödinger equation compatible with special relativity). One can note that the ordinary Schrödinger (free) equation
$$\hat{H}\psi=ih\frac{\partial \psi}{\partial t}$$ 
with a hamiltonian of the form $H=\hat{p}^2/2m$ reproduces the classical energy-momentum relation $E=p^2/2m$ which is not relativistic.
One then can use the same association of canonical variables into operators but in the relativistic energy-momentum relation $E^2=m^2 + p^2$ (1) (taking $c=1$), in order to derive a relativistic version of Schrödinger equation.
If one does that in relation (1), one gets the Klein-Gordon equation. This equation in good for bosons and the energy is bounded from below when bosons are described by this equation.
However if one tries to apply the KG equation to fermions, problems arise. For starters, fermions which obey Pauli's principle (thus have to be quantised with anti-commutators), make the KG equation have unbounded energy from below (which is something equivalent to a perpetuum mobile).
The problem is that the relation (1) is second order in the Energy (thus in equation form, is second order in time evolution $\partial/\partial t$).
Dirac tried to solve this by factoring the relation (1) into this form (2):
$$E=\mathbf{\alpha} m + \mathbf{\beta} p$$
but in order for the relation (2) to be equal to relation (1) when squared, $\mathbf{\alpha}$ and $\mathbf{\beta}$ are not ordinary numbers but some kind of matrices (specificaly pauli matrices). Thus we enter what is refered as formalism of spinors in the Dirac equation which is a quantum equation derived from relation (2).
Indeed the Dirac equation can describe fermions (aka Pauli exclusion principle, anti-commutators) and the energy is in fact bounded from below.
One artifact is that the dirac equation describes 2 particles and not one, since matrices/spinors are involved. Eventualy this led to the discovery of the positron (and anti-particles/anti-matter). 
Finaly, the synthesis of quantum mechanics and relativity (along with anti-particles) led to the Quantum Field Theory formalism and the Spin-Statistics theorem, which relates theoreticaly the spin of a particle with the type of statistics it follows (integral spin->Bose-Einstein statistics, half-integral spin->Fermi-Dirac statistics)
References:


*

*PAM Dirac, The Quantum Theory of the Electron

*PAM Dirac, The Quantum Theory of the Electron. Part II

*PAM Dirac, A Positive-Energy Relativistic Wave Equation

*PAM Dirac, A Positive-Energy Relativistic Wave Equation. II
A: Spin-1/2 admits first order equations simply because 
$
(\mathbf{1/2,1/2})\otimes (\mathbf{0,1/2})
$
contains the representation $(\mathbf{1/2,0})$ so that a linear equation for free particles can be written (i.e. it contains a derivative acting on one field and returning one field).
The first term in the product is the derivative that transforms as a Lorentz vector $(\mathbf{1/2,1/2})$, whereas $(\mathbf{0,1/2})$ and $(\mathbf{1/2,0})$ are left and right-handed spinors respectively. The clebsch-Gordan coefficients are nothing but that the gamma matrices. 
Clearly the same is forbidden for the integer spins. For a scalar is trivial. For a spin-1 field $(\mathbf{1/2,1/2})$ one has that $( \mathbf{ 1/2, 1/2}) \otimes ( \mathbf{1/2,1/2}) $ does not contain $(\mathbf{1/2,1/2})$, and so on for higher integer spins. It is basically the group structure of Lorentz symmetry that forbids first order eq. for integer spins.
A: Fermions obey the Fermi-Dirac statistics, while Bosons the Bose-Einstein statistics. This is an experimental fact and we can't do anything about it. You find its first and most famous evidence in the Pauli principle. To mention some more, Bose condensation and Fermi blocking are a fact of everyday science, and we have even confirmation that the wavefunction of a fermion changes its sign after a 360 degree rotation around. These things experimentally distinguish fermions from bosons. 
The Dirac equation brings all the fermionic features in the game. It does it in a first order differential equation. On the other hand, the Klein Gordon equation brings all the bosonic features in the game. It does this in second order. You cannot get fermionic features from the Klein Gordon equation or vice verse. On this, there is a nice paragraph in Peskin&Schroeder book at chapter 3.5: "How not to quantize the Dirac field: a lesson in Spin and Statistics." It shows that something goes terribly wrong if you try to quantize the Dirac equation in the same manner you quantize bosonic particles. 
Statistics is what we observe, statistics is what we try to model. 
A: Steven Weinberg starts with symmetries of relativity, and quantum framework, and arrives in chapter 5 of volume one at Dirac equation. He first gets the quantum field for spin 1/2 particles (without reference to any Lagrangian or wave equation), then he constructs vacuum expectation value of the field and its 
adjoint, this leads to Feynman-Dyson propagator. Earlier in the analysis he gets the Dirac operator when he calculates spin sums. Acting the Dirac operator on the FD propagator then yields, up to a sign, the four dimensional Dirac delta function for x-y.  From this follows the Lagrangian density.
Mathematically, the linearity depends on the fact that (\vec Pauli matrix . unit momentum vector)^2 = identity matrix. 
If you need more details, please consult my monograph (Mass Dimension One Fermions, Cambridge Monographs on Mathematical Physics, Cambridge University Press, 2019). And one may evade the linearity you write about.
Wait for a preprint soon -- on spin 1/2 bosons with linear wave equation.
A: In fact, the Dirac equation can be expressed as a second-order differential equation of the form
\begin{equation}
\frac{1}{c^2}\frac{\partial^2\psi}{\partial t^2}-\nabla^2\psi\pm\frac{2mi}{\hbar}\frac{\partial\psi}{\partial t}-\left(\frac{mc}{\hbar}\right)^2\psi=0\,.
\end{equation}
