We study the effects of anisotropic thermal conduction on magnetized spherical accretion flows using global axisymmetric MHD simulations. In low collisionality plasmas, the Bondi spherical accretion solution is unstable to the magnetothermal instability (MTI). The MTI grows rapidly at large radii where the inflow is subsonic. For a weak initial field, the MTI saturates by creating a primarily radial magnetic field, i.e., by aligning the field lines with the background temperature gradient. The saturation is quasilinear in the sense that the magnetic field is amplified by a factor of 10-30 independent of the initial field strength (for weak fields). In the saturated state, the conductive heat flux is much larger than the convective heat flux, and is comparable to the field-free (Spitzer) value (since the field lines are largely radial). The MTI by itself does not appreciably change the accretion rate \dot M relative to the Bondi rate \dot MB. However, the radial field lines created by the MTI are amplified by flux freezing as the plasma flows in to small radii. Oppositely directed field lines are brought together by the converging inflow, leading to significant resistive heating. When the magnetic energy density is comparable to the gravitational potential energy density, the plasma is heated to roughly the virial temperature; the mean inflow is highly subsonic; most of the energy released by accretion is transported to large radii by thermal conduction; and the accretion rate \dot M \ll \dot MB. The predominantly radial magnetic field created by the MTI at large radii in spherical accretion flows may account for the stable Faraday rotation measure towards Sgr A* in the Galactic Center.