Let $(M, g)$ be a complete Riemannian manifold, $L=\\Delta -\\nabla \\phi \\cdot \\nabla$ be a Markovian symmetric diffusion operator with an invariant measure $d\\mu(x)=e^{-\\phi(x)}d\\nu(x)$, where $\\phi\\in C^2(M)$, $\\nu$ is the Riemannian volume measure on $(M, g)$. A fundamental question in harmonic analysis and potential theory asks whether or not the Riesz transform $Ra(L)=\\nabla(a-L)^{-1/2}$ is bounded in $L^p(\\mu)$ for all $1<p<\\infty$ and for certain $a\\geq 0$. An affirmative answer to this problem has many important applications in elliptic or parabolic PDEs, potential theory, probability theory, the $L^p$-Hodge decomposition theory and in the study of Navier-Stokes equations and boundary value problems. Using some new interplays between harmonic analysis, differential geometry and probability theory, we prove that the Riesz transform $R_a(L)=\\nabla(a-L)^{-1/2}$ is bounded in $L^p(\\mu)$ for all $a>0$ and $p\\geq 2$ provided that $L$ generates a ultracontractive Markovian semigroup $P_t=e^{tL}$ in the sense that $P_t 1=1$ for all $t\\geq 0$, $|P_t|{1, \\infty} < Ct^{-n/2}$ for all $t\\in (0, 1]$ for some constants $C>0$ and $n > 1$, and satisfies $$ (K+c)^{-}\\in L^{{n\\over 2}+\\epsilon}(M, \\mu) $$ for some constants $c\\geq 0$ and $\\epsilon>0$, where $K(x)$ denotes the lowest eigenvalue of the Bakry-Emery Ricci curvature $Ric(L)=Ric+\\nabla^2\\phi$ on $T_x M$, i.e., $$ K(x)=\\inf\\limits{Ric(L)(v, v): v\\in T_x M, |v|=1}, \\quad\\forall\\ x\\in M. $$ Examples of diffusion operators on complete non-compact Riemannian manifolds with unbounded negative Ricci curvature or Bakry-Emery Ricci curvature are given for which the Riesz transform $R_a(L)$ is bounded in $L^p(\\mu)$ for all $p\\geq 2$ and for all $a>0$ (or even for all $a\\geq 0$).