Let $J$ be the Lévy density of a symmetric Lévy process in $\mathbb{R}^d$ with its Lévy exponent satisfying a weak lower scaling condition at infinity. Consider the non-symmetric and non-local operator $$ {\mathcal L}^κf(x):= \lim_{ε\downarrow 0} \int_{\{z \in \mathbb{R}^d: |z|>ε\}}(f(x+z)-f(x))κ(x,z)J(z)\, dz\, , $$ where $κ(x,z)$ is a Borel measurable function on $\mathbb{R}^d\times \mathbb{R}^d$ satisfying $0<κ_0\le κ(x,z)\le κ_1$, $κ(x,z)=κ(x,-z)$ and $|κ(x,z)-κ(y,z)|\le κ_2|x-y|^β$ for some $β\in (0, 1)$. We construct the heat kernel $p^κ(t, x, y)$ of ${\mathcal L}^κ$, establish its upper bound as well as its fractional derivative and gradient estimates. Under an additional weak upper scaling condition at infinity, we also establish a lower bound for the heat kernel $p^κ$.