Why propagator is different with different spin? I could not understand what does propagators do with the spin.
Is it simply because the Lagrangian is different with different spins?
 A: This is easy to understand from the functional approach. Let me outline it here . First a small disclaimer: there might be some incorrect numeric factors, so if they are spotted let me know. Consider a free theory, which therefore has a quadratic action $S_0[\Phi]$ where $\Phi$ can be any field. Let us suppose that $\Phi$ is real and takes values in some vector space $V$, so that we write its components $\Phi^I(x)$ in some basis. The case for complex fields is completely analogous and I'll ommit it.
Now, it is always possible to write $S_0[\Phi]$ as $$S_0[\Phi]=-\dfrac{1}{2}\int d^Dx d^Dy\ \Phi^I(x){\cal K}^{I}_{\phantom IJ}(x,y)\Phi^J(y)\tag{1}$$
where ${\cal K}^{I}_{\phantom IJ}(x,y)$ is one differential operator viewed here as a matrix in field space. It maps $$\Phi^I(x)\to ({\cal K}\Phi)^I(x)=\int d^Dy {\cal K}^{I}_{\phantom I J}(x,y)\Phi^J(y)\tag{2}.$$
For example, if we have a scalar field $\phi$ then ${\cal K}=\Box+m^2$ is just the KG operator. In general ${\cal K}$ is the operator such that the free wave equation that $\Phi$ obeys is ${\cal K}\Phi=0$.
That said, let $\Delta^I_{\phantom{I}J}(x,y)$ be an inverse of ${\cal K}^{I}_{\phantom IJ}(x,y)$ defined by the Green's function differential equation $$\int d^Dy{\cal K}^I_{\phantom IJ}(x,y)\Delta^J_{\phantom J K}(y,z)=\delta^I_{\phantom IK}\delta^{(D)}(x-z)\tag{3}.$$
In Lorentzian signature (3) must be supplemented by further conditions to make the inverse $\Delta^I_{\phantom{I}J}(x,y)$ well-defined: this is just the $i\epsilon$ prescription, but I'll get to it. The basic point to observe is that the problem is that there is no unique choice for such an inverse. So for the moment assume we have picked one.
Now the basic point of all of this is that the study of Gaussian functional integrals tells us that the generating functional of correlation functions associated to (1) is given by $$Z_0[J]=\int \mathfrak{D}\Phi e^{iS_0[\Phi]+i\int d^Dx J_I(x)\Phi^I(x)}:=\exp\left[\frac{i}{2}\int d^Dx d^Dy 
 J^I(x)\Delta_{I J}(x,y)J^J(y)\right]\tag{4}.$$
Once $Z_0[J]$ is known the correlation functions follow from functional differentiation. In particular, the two-point function is $$\langle \Phi^I(x)\Phi^J(y)\rangle:=(-i)^2\dfrac{\delta^2 Z_0[J]}{\delta J_I(x)\delta J_J(y)}\bigg|_{J=0}=\Delta^{IJ}(x,y)\tag{5}.$$
So the inverse we have chosen will be the two-point function, which is exactly what we define as the propagator. If we want to generate time-ordered two-point functions this imposes a condition on the solution of (3) which turns into the $i\varepsilon$ prescription and uniquely defines $\Delta^I_{\phantom I J}(x,y)$.
The same story, mutatis mutandis, can be developed for the complex and fermionic case.
So how this answers your question? Basically, the propagator follows directly as an inverse of the wave operator appearing in the classical action. Such an operator depends on spin. For spin zero we have the KG operator, for spin half we have the Dirac (or Weyl) operator, and so forth. Since the operator changes for the various fields, the propagator changes as well.
