  In $\mathbb R^d$, $d \geq 3$, consider the divergence and the non-divergence
form operators \begin{equation} \tag{$i$} - \nabla \cdot a \cdot \nabla + b
\cdot \nabla, \end{equation} \begin{equation} \tag{$ii$} - a \cdot \nabla^2 + b
\cdot \nabla, \end{equation} where $a=I+c \mathsf{f} \otimes \mathsf{f}$, the
vector fields $\nabla_i \mathsf{f}$ ($i=1,2,\dots,d$) and $b$ are form-bounded
(this includes the sub-critical class $[L^d + L^\infty]^d$ as well as vector
fields having critical-order singularities). We characterize quantitative
dependence on $c$ and the values of the form-bounds of the $L^q \rightarrow
W^{1,qd/(d-2)}$ regularity of the resolvents of the operator realizations of
($i$), ($ii$) in $L^q$, $q \geq 2 \vee ( d-2)$ as (minus) generators of
positivity preserving $L^\infty$ contraction $C_0$ semigroups. The latter
allows to run an iteration procedure $L^p \rightarrow L^\infty$ that yields
associated with ($i$), ($ii$) $L^q$-strong Feller semigroups.
