The goal of this paper is to establish sharp two-sided estimates on the heat kernels of two types of purely discontinuous symmetric Markov processes in the upper half-space of $\mathbb R^d$ with jump kernels degenerate at the boundary. The jump kernels are of the form $J(x,y)=\mathcal B(x,y)|x-y|^{-α-d}$, $α\in (0,2)$, where the function $\mathcal B$ depends on four parameters and may vanish at the boundary. Our results are the first sharp two-sided estimates for the heat kernels of non-local operators with jump kernels degenerate at the boundary. The first type of processes are conservative Markov processes on $\overline{\mathbb R}^d_+$ with jump kernel $J(x,y)$. Depending on the regions where the parameters belong, the heat kernels estimates have three different forms, two of them are qualitatively different from all previously known heat kernel estimates. The second type of processes are the processes above killed either by a critical potential or upon hitting the boundary of the half-space. We establish that their heat kernel estimates have the approximate factorization property with survival probabilities decaying as a power of the distance to the boundary, where the power depends on the constant in the critical potential.