In this paper, a high-order solver combining the flux reconstruction (FR) method and the thermal lattice Boltzmann flux solver (FRTLBFS) is developed for accurately and efficiently simulating incompressible thermal flows. The conservative differential equations recovered from Chapman-Enskog analysis of the thermal lattice Boltzmann equation are solved by the high-order FR method. The thermal lattice Boltzmann method is only applied to reconstruct the local solution used for evaluating fluxes at the solution and flux points. Unlike the traditional Navier-Stokes-Boussinesq (NSB) solvers where the inviscid and viscous terms are treated separately, the inviscid and viscous fluxes in the current FRTLBFS are coupled and computed uniformly. In comparison with the recently developed high-order flux reconstruction thermal lattice Boltzmann method, the FRTLBFS holds advantages such as high-order accuracy, good stability, and compactness but is more efficient and low storage, since only macroscopic flow variables including density, velocity, and temperature are stored and evolved. In addition, the physical boundary conditions in FRTLBFS can be directly implemented by using the same method as in conventional NSB solvers. Numerical validations of the proposed method are implemented by simulating (a) the porous plate problem, (b) natural convection in a square cavity, (c) unsteady natural convection in a tall cavity, and (d) thermal lid-driven cavity flow. Numerical results demonstrate that the present solver is an attractive tool to simulate incompressible thermal flows due to its high-order accuracy, stability, and low memory cost.