Let $d\ge1$ and $0<α<2$. Consider the integro-differential operator \[ \mathcal{L}f(x) =\int_{\mathbb{R}^{d}\backslash\{0\}}\left[f(x+h)-f(x)-χ_α(h)\nabla f(x)\cdot h\right]\frac{n(x,h)}{|h|^{d+α}}\mathrm{d}h+\mathbf{1}_{α>1}b(x)\cdot\nabla f(x), \] where $χ_α(h):=\mathbf{1}_{α>1}+\mathbf{1}_{α=1}\mathbf{1}_{\{|h|\le1\}}$, $b:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is bounded measurable, and $n:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ is measurable and bounded above and below respectively by two positive constants. Further, we assume that $n(x,h)$ is Hölder continuous in $x$, uniformly with respect to $h\in\mathbb{R}^{d}$. In the case $α=1,$ we assume additionally $\int_{\partial B_{r}}n(x,h)h\mathrm{d}S_{r}(h)=0$, $\forall r \in (0,\infty)$, where $\mathrm{d}S_{r}$ is the surface measure on $\partial B_{r}$, the boundary of the ball with radius $r$ and center $0$. In this paper, we establish two-sided estimates for the heat kernel of the Markov process associated with the operator $\mathcal{L}$. This extends a recent result of Z.-Q. Chen and X. Zhang.