A Spectral Difference (SD) algorithm on tensor-product elements which solves the reacting compressible Navier-Stokes equations is presented. The classical SD algorithm is shown to be unstable when a multi-species gas where thermodynamic properties depend on temperature and species mass fractions is considered. It is demonstrated that computing primitive variables from conservative variables at solution points and extrapolating them at flux points, rather than conservative variables, make the SD method stable for such case. In addition, a new methodology, more stable, for computing the diffusive fluxes at an interface is detailed. To allow the simulation of combustion, characteristic and wall boundary conditions for multi-species flows in the SD framework are introduced and the thickened flame turbulent combustion model is adapted to the SD context. Validation test cases from one-dimensional and two-dimensional Flow, Turbulence and Combustion laminar flames to a three-dimensional turbulent flame have been performed showing excellent agreement with each of the reference solutions. To the authors' knowledge, this is the first time that a 3D turbulent flame is simulated using the SD method. The interest of employing a high-order method, such as the SD, in combustion is demonstrated: high values of the polynomial degree improve the quality of the results which advocates for the use of p-refinement when simulating reacting flows.