  In this paper, we extend the operator-split asymptotic-preserving,
semi-Lagrangian algorithm for time dependent anisotropic heat transport
equation proposed in [Chac\'on et al., JCP, 272, 719-746, 2014] to use a fully
implicit time integration with backward differentiation formulas. The proposed
implicit method can deal with arbitrary heat-transport anisotropy ratios
\chi_{\parallel}/\chi_{\perp}\ggg1 (with \chi_{\parallel}, \chi_{\perp} the
parallel and perpendicular heat diffusivities, respectively) in complicated
magnetic field topologies in an accurate and efficient manner. The implicit
algorithm is second-order accurate temporally, has favorable positivity
preservation properties, and demonstrates an accurate treatment at boundary
layers (e.g., island separatrices), which was not ensured by the operator-split
implementation. The condition number of the resulting algebraic system is
independent of the anisotropy ratio, and is inverted with preconditioned GMRES.
We propose a simple preconditioner that renders the linear system compact,
resulting in mesh-independent convergence rates for topologically simple
magnetic fields, and convergence rates scaling as \sim(N\Delta t)^{1/4} (with N
the total mesh size and \Delta t the timestep) in topologically complex
magnetic-field configurations. We demonstrate the accuracy and performance of
the approach with test problems of varying complexity, including an
analytically tractable boundary-layer problem in a straight magnetic field, and
a topologically complex magnetic field featuring magnetic islands with extreme
anisotropy ratios (\chi_{\parallel}/\chi_{\perp}=10^{10}).
