Let $L:= -a(x) (-Δ)^{α/2}+ (b(x), \nabla)$, where $α\in (0,2)$, and $a:\rd\to (0,\infty)$, $b: \rd\to \rd$. Under certain regularity assumptions on the coefficients $a$ and $b$, we associate with the $C_\infty(\rd)$-closure of $(L, C_\infty^2(\rd))$ a Feller Markov process $X$, which possesses a transition probability density $p_t(x,y)$. To construct this transition probability density and to obtain the two-sided estimates on it, we develop a new version of the parametrix method, which allows us to handle the case $0<α\leq 1$ and $b\neq 0$, i.e. when the gradient part of the generator is not dominated by the jump part..