In this paper, we consider the following type of non-local
(pseudo-differential) operators $\mathcal{L}$ on $\mathbb{R}^d$:
\begin{align*}
\mathcal{L} u(x) = \frac{1}{2} \sum_{i, j=1}^d \frac{\partial}{\partial x_i}
\Big(a_{ij}(x) \frac{\partial u(x)}{\partial x_j}\Big)
\\+
\lim_{\varepsilon \downarrow 0} \int_{\{y\in \mathbb{R}^d: |y-x|>\varepsilon\}}
(u(y)-u(x)) J(x, y) dy,
\end{align*}
where $A(x)=(a_{ij}(x))_{1\leq i,j\leq d}$ is a measurable $d\times d$ matrix-valued
function on $\mathbb{R}^d$ that is uniformly
elliptic and bounded and $J$ is a symmetric measurable non-trivial
non-negative kernel on $\mathbb{R}^d \times \mathbb{R}^d$ satisfying certain
conditions. Corresponding to $\mathcal{L}$ is a symmetric strong Markov
process $X$ on $\mathbb{R}^d$ that has both the diffusion component and pure
jump component. We establish a priori Hölder estimate for bounded
parabolic functions of $\mathcal{L}$ and parabolic Harnack principle for
positive parabolic functions of $\mathcal{L}$. Moreover, two-sided sharp
heat kernel estimates are derived for such operator $\mathcal{L}$ and
jump-diffusion $X$. In particular, our results apply to the mixture
of symmetric diffusion of uniformly elliptic divergence form
operator and mixed stable-like processes on $\mathbb{R}^d$. To establish
these results, we employ methods from both probability theory and
analysis.