  The purpose of this paper is to introduce and study Poincar\'e-Steklov (PS)
operators associated to the Dirac operator $D_m$ with the so-called MIT bag
boundary condition. In a domain $\Omega\subset\mathbb{R}^3$, for a complex
number $z$ and for $U_z$ a solution of $(D_m-z)U_z=0$, the associated PS
operator maps the value of $\Gamma_- U_z$, the MIT bag boundary value of $U_z$,
to $\Gamma_+ U_z$, where $\Gamma_\pm$ are projections along the boundary
$\partial\Omega$ and $(\Gamma_ - + \Gamma_+) = t_{\partial\Omega}$ is the trace
operator on $\partial\Omega$.
  In the first part of this paper, we show that the PS operator is a zero-order
pseudodifferential operator and give its principal symbol. In the second part,
we study the PS operator when the mass $m$ is large, and we prove that it fits
into the framework of $1/m$-pseudodifferential operators, and we derive some
important properties, especially its semiclassical principal symbol.
Subsequently, we apply these results to establish a Krein-type resolvent
formula for the Dirac operator $H_M= D_m+ M\beta
1_{\mathbb{R}^3\setminus\overline{\Omega}}$ for large masses $M>0$, in terms of
the resolvent of the MIT bag operator on $\Omega$. With its help, the large
coupling convergence with a convergence rate of $\mathcal{O}(M^{-1})$ is shown.
