AbstractLet {L=-\operatorname{div}(A\nabla)} be a second-order divergence form elliptic operator and let A be an accretive, {n\times n} matrix with bounded measurable complex coefficients in {{\mathbb{R}}^{n}}. Let {L^{-\frac{\alpha}{2}}} be the fractional integral associated to L for {0<\alpha<n}. For {b\in L_{\mathrm{loc}}({\mathbb{R}}^{n})} and {k\in{\mathbb{N}}}, the k-th order commutator of b and {L^{-\frac{\alpha}{2}}} is given by(L^{-\frac{\alpha}{2}})_{b,k}f(x)=L^{-\frac{\alpha}{2}}((b(x)-b)^{k}f)(x).In the paper,
we mainly show that if {b\in\mathrm{BMO}({\mathbb{R}}^{n})}, {0<\lambda<n} and {0<\alpha<n-\lambda}, then {(L^{-\frac{\alpha}{2}})_{b,k}} is bounded from {L^{p,\lambda}} to {L^{q,\lambda}} for {p_{-}(L)<p<q<p_{+}(L)\frac{n-\lambda}{n}} and {\frac{1}{q}=\frac{1}{p}-\frac{\alpha}{n-\lambda}}, where {p_{-}(L)} and {p_{+}(L)} are the two critical exponents for the {L^{p}} uniform
boundedness of the semigroup {\{e^{-tL}\}_{t>0}}. Also, we establish the boundedness of the commutator of the fractional integral with Lipschitz function on Morrey spaces. The results encompass
what is known for the classical Riesz potentials and elliptic operators with Gaussian domination by the classical heat
operator.