Dirac wanted to reformulate a relativistic version of the Schrödinger equation without the problems arising in the Klein-Gordon equation (for example, the fact that the probability density is negative in the KG equation). So he modified the KG equation
$$(\partial_\mu\partial^{\mu} + m^2)\psi=0$$
to the first order equation
$$(i\gamma^\mu\partial_{\mu}-m)\psi=0$$
which can be proved that gives a positive definite probability density. But what are these coefficients $\gamma^\mu$? Applying $i\gamma^{\nu}\partial_{\nu}$ on both sides of the equation gives
$$-\gamma^{\nu}\gamma^{\mu}\partial_{\nu}\partial_{\mu}\psi=im\gamma^{\nu}\partial_{\nu}\psi=m^2\psi$$
Where in the last step we have used the result of the original Dirac equation. Now, we can write the last equation as:
$$-\frac{1}{2}(\gamma^{\nu}\gamma^{\mu}+\gamma^{\mu}\gamma^{\nu})\partial_{\nu}\partial_{\mu}\psi=m^2\psi$$
We see that the term in the parenthesis is the well known anticommutator $\{ \gamma^{\nu},\gamma^{\mu} \}=\gamma^{\nu}\gamma^{\mu}+\gamma^{\mu}\gamma^{\nu}$. We know as well that if we have a plain wave $\psi\sim e^{ip^{\mu}x_{\mu}}$ and since $p_{\mu}p^{\mu}=m^2$ the Dirac equation will become:
$$\frac{1}{2}\{ \gamma^{\nu},\gamma^{\mu} \}p_{\nu}p_{\mu}=m^2$$
It's easy to see that in order to fulfill the previous equations, the anticommutator has to be:
$$\{ \gamma^{\nu},\gamma^{\mu} \}=2\eta^{\nu\mu}$$
Since $\eta^{\nu\mu}$ are the components of a $4\times4$ matrix, we can easily deduce that the four quantities $\gamma^\mu$ have to be matrices as well. The above relation already defines the $\gamma^{\mu}$, but we will see some other nice properties. The easy ones are:
$$(\gamma^0)^2=\eta^{00}=1, \ \ (\gamma^i)^2=-1, \ \ \gamma^0\gamma^i+\gamma^i\gamma^0=0, \ \ \gamma^i\gamma^j+\gamma^j\gamma^i=2\delta^{ij}$$
From the first one we can deduce that the eigenvalues of $\gamma^0$ are $+1$ and $-1$. We have also that if $\nu\neq\mu$:
$$\gamma^{\mu}\gamma^{\nu}=-\gamma^{\nu}\gamma^{\mu}$$
In particular, if we take $\nu=0$ and $\mu=i$ we get:
$$\gamma^0\gamma^i=-\gamma^i\gamma^0, \ \text{multiplying by $\gamma^i$ on the right}\to\gamma^0=\gamma^i\gamma^0\gamma^i$$
If we take the trace of the last result:
$$\text{Tr}(\gamma^0)=\text{Tr}(\gamma^i\gamma^0\gamma^i)=\text{Tr}(\gamma^i\gamma^i\gamma^0)=\text{Tr}(-\gamma^0)=-\text{Tr}(\gamma^0)\Longrightarrow \text{Tr}(\gamma^0)=0$$
Here we have used some properties of the trace, namely $\text{Tr}(AB)=\text{Tr}(BA)$ and $\text{Tr}(-A)=-\text{Tr}(A)$. From this fact and that the eigenvalues of $\gamma^0$ are $+1$ and $-1$ we can see that this matrix has to be of even dimension. It can't be of dimension 2 because it would not fulfill the above properties. It's not hard to see that its dimension has to be 4. Now we can look upon the Dirac equation:
$$i\gamma^{\mu}\partial_{\mu}\psi=m\psi\Longrightarrow i\gamma^{0}\frac{\partial\psi}{\partial t}=i\gamma^i\frac{\partial\psi}{\partial x^i}+m\psi$$
This equation only makes sense if $\psi$ is a 4-dimensional object, because we have $4\times 4$ matrices acting on it. We call $\psi$ a spinor (4-spinor, Dirac spinor), and we will see that is defined by its transformation properties under Lorentz transformations. Recalling the Dirac equation and multiplying by $\gamma^0$ on both sides:
$$i\frac{\partial\psi}{\partial t}=i\gamma^0\gamma^i\frac{\partial\psi}{\partial x^i}+m\gamma^0\psi$$
We can identify the left hand side of the equation as the hamiltonian operator $H\psi=i\frac{\partial\psi}{\partial t}$, and obviously has the form:
$$H=\gamma^0\gamma^ip_i+m\gamma^0$$
Sometimes it's convenient to express it as:
$$ H=\boldsymbol{\alpha}\cdot\mathbf{p}+\beta m$$
Where we define the operators $\boldsymbol{\alpha}$ and $\beta$:
$$\boldsymbol{\alpha}=\left(\begin{array}{cc}
0 & \boldsymbol{\sigma} \\
\boldsymbol{\sigma} & 0
\end{array}\right), \ \ \beta=\gamma^{0}$$
By requiring to be hermitian $H=H^{\dagger}$ we have:
$$(\gamma^0\gamma^ip_i)^{\dagger}+(m\gamma^0)^{\dagger}=\gamma^0\gamma^ip_i+m\gamma^0$$
Since the momentum is already hermitean $p_i=p_i^{\dagger}$ we have the conditions:
$$(\gamma^0)^{\dagger}=\gamma^0, \ \ (\gamma^0\gamma^i)^{\dagger}=\gamma^0\gamma^i$$
From the first one and since $(\gamma^0\gamma^i)^{\dagger}=(\gamma^i)^{\dagger}(\gamma^0)^{\dagger}$ we see that $(\gamma^i)^{\dagger}=-\gamma^i$. One can write all this properties in a single formula:
$$\gamma^0(\gamma^\mu)^{\dagger}\gamma^0=\gamma^{\mu}$$
Indeed there are not only one set of matrices that fulfill the above properties. A common way of writing them is what is known as the Dirac representation:
$$
\gamma^0=\left(\begin{array}{cc}
\mathbb{I} & 0 \\
0 & - \mathbb{I}
\end{array}\right), \ \
\gamma^i=\left(\begin{array}{cc}
0 & \sigma^i \\
-\sigma^i & 0
\end{array}\right)
$$
Where of course $\sigma^i$ are the Pauli matrices:
$$\sigma^1=\left(\begin{array}{cc}
0 & 1 \\
1 & 0
\end{array}\right), \
\sigma^2=\left(\begin{array}{cc}
0 & -i \\
i & 0
\end{array}\right), \
\sigma^3=\left(\begin{array}{cc}
1 & 0 \\
0 & -1
\end{array}\right)$$
Therefore, the Pauli matrices appear as a result of several conditions the $\gamma$ matrices must fulfill. Some other well known representations of the $\gamma$ matrices that fulfill the properties are the Weyl representation or the Majorana representation, both equally useful in their respective aspects. In fact, all the representations are equivalent up to a unitary transformation:
$$\gamma^{\mu}_{_{R'}}=\mathcal{U}\gamma^{\mu}_{_R}\mathcal{ \ U}^{-1}$$
Where $\mathcal{U}$ is a unitary operator and $R$ and $R'$ label the representation.
Sorry for the long explanation, but I hope it helps!
